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ABSTRACT 


The accumulation of static charge on the surface of spacecraft 
can produce unknown potential fields that may cause error in the 
measurement of scientific data. Therefore methods for determining 
potential fields are being sought. This report describes a general 
technique for expanding an unknown potential field in terms of a linear 
summation of weighted dipole or quadrupole fields. Laplace's equation 
describes the unknown potential in regions that are devoid of electric 
charge when the nearby surface charge distribution is assumed to be 
independent of time. Therefore, general solutions to Laplace's equa- 
tion can be generated by summing multipole fields (i.e., particular 
solutions to Laplace's equation) as long as the multipole locations 
where the multipole fields have singularities are not included in the 
spatial region of interest. The general solution approximating an un- 
known potential near a charged surface can be developed if suitable 
boundary conditions are available or if measurements of charged 
particle trajectories in the unknown potential can be obtained. 

The classical boundary value problem in electrostatics 
can be described as finding the solution to Laplace's equation when 
the electric potential is known on a boundary enclosing a spatial do- 
main of interest. In this report, computational methods are developed 
for the iterative addition of dipole fields until solutions to the 
classical boundary value problem can be obtained. Various solution 
potentials are compared inside the boundary with a more precise cal- 
culation of the potential to derive optimal schemes for locating the 
singularities of the dipole fields (e.g., dipoles should not be placed 
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too close to the boundary) . Then, the problem of determining solutions 
to Laplace' s equation on an unbounded domain as constrained by pertinent 
electron trajectory data is considered. The initial electron coordinates 
and velocities, as well as the final electron coordinates, comprise a 
set of constraints on the various schemes that are developed for dipole 
and quadrupole synthesis of approximations to the test potential. The 
various schemes are then compared in terms of convergence limits and 
rates and in terms of their accuracy on a finite domain. The report is 
concluded with a description of an electron-beam apparatus used for 
making trajectory measurements in a bell jar. The best computational 
schemes are used to synthesize approximations of a simulated test poten- 
tial created in the bell jar. 
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I. INTRODUCTION 


The accumulation of electrostatic charge on the surface of a 

spacecraft produces an electric potential field in nearby spatial re- 

2 

gions. Laplace’s equation, V ip = 0, defines a class of functions which 

may represent this potential field everywhere outside the spacecraft 

(i.e., spatial regions having little charge). If the surface charge is 

known, the potential field ip can be determined by an integration over 

the spacecraft of the field due to a point charge weighted with the 

known surface charge distribution. However, for many applications the 

surface charge is unknown and must be determined by solving Laplace's 

equation constrained by the pertinent boundary conditions that may be 

known. The potential fields of the ideal dipole, quadrupole, octapole, 

etc. are solutions of Laplace’s equations in regions bounded away from 

the sites of these ideal poles. Therefore, the unknown potential ip 

can be synthesized as a linear sum of multipole potential fields, with 

all point multipoles located outside the region of interest (ref. 1). 

The classical boundary value problem may be considered for 

Laplace’s equation in terms of a potential (}> unknown in region R, but 

known at a finite number of points on a simple closed curve containing 

the region R. In two dimensions, the method of finite differences 

could be employed to solve for the unknown potential at the nodes of 

a grid covering region R subject to the known boundary values of the 

potential ij; (e.g., the specification of the potential at 4N points on 

2 

a square boundary determines an approximation of the potential at N 
points on the interior of the square boundary). Consider instead, 
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an iterative solution technique. Suppose some potential (j) exists as an 
estimate of ip. Then a measure of the discrepancy D between ip and (p 
can be defined as the sxim of the absolute differences between \p and (j) 
at the finite number of nodes on the boundary where ij; is known. The 
estimated potential (|) can be said to be converging toward ijj if the dis- 
crepancy D is reduced upon successive applications of the iterative 
scheme which may include operations of rotation, shifting or scaling. 
However, discrepancy reduction and the shaping of the potential approxi- 
mation will be primarily accomplished by the summation of multipole 
fields successively centered at points outside the boundary near loca- 
tions where the discrepancy D is largest. 

The expansion of potential fields is most often expressed as the 
linear sum of weighted multipole fields (e.g., one each of monopole, 
dipole, quadrupole, etc.) with all pole locations superimposed. This 
practice produces only one singularity in the spatial solution space. 
However, we have adopted the practice of allowing multipole locations 
along a line parallel to the surface and at a specified distance from 
it (i.e., usually all alike dipoles or quadrupoles) . Thus, field shaping 
is facilitated but many field singularities are produced. 

When considering the potential field existing near a charged sur- 
face, the potential may be unknown in an unbounded spatial region. 

Due to insufficient boundary conditions, additional constraints must be 
specified in order to approximate the potential on a finite spatial re- 
gion. Actually, in mathematical terms, the region is bounded at 
infinity where the potential must equal zero. In order to obtain 
measurements of potential fields existing near charged surfaces. 
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consider an electron beam emanating from a source point (X0,Y0) on a 
zero potential reference plane. If the- reference plane were 
located far from the charged surface, it would approximate the boundary 
condition at infinity. Practically speaking, the reference plane can 
be moved to a finite and measurable distance from the charged surface 
while producing only minimal distortions of the potential field. If 
the spacecraft surface is negatively charged the electron beam will be 
deflected back toward the reference plane where the coordinates of a 
sink point (XT,YT) can be recorded. As seen in Figure 1, a displace- 
ment vector measured from the beam source point to the beam sink point 
represents one constraint on the family of functions that satisfy 
Laplace's equation. Suppose that a potential field, partially bounded 
by a charged surface and the reference plane, is to be determined based 
on measurements of many such displacement vectors associated with 
electron beams with known energies and initial directions. Suppose 
further that some estimate <() of the potential if; exists on the region 
of interest. Estimated displacement vectors can be obtained from the 
initial electron beam data and the estimated potential by programming 
the laws of electron mechanics. These estimated displacements can be 
compared to those measured as data to define the error related to the 
electron trajectory constraints. A discrepancy DD can be defined as 
the sum of the absolute differences between the estimated displacement 
vectors and the measured displacement vectors. The estimated potential 
(t) is said to be converging toward i(; if the discrepancy is reduced upon 
successive applications of some iterative scheme which may include 
scaling or rotating operations. However, major alterations in the 
estimated potential can be accomplished by adding the potential fields 
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Figure 1. The Orientation of Ideal Line 
Multipoles with Respect to a Charged 
Surface and Reference Plane 
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of pairs of ideal multipoles. One multipole can be located above the 
charged surface at a given distance from the reference plane, while 
another multipole is located an equal distance below the reference 
plane. This imaging technique maintains a potential of zero on the 
reference plane while allowing for the shaping of the estimated poten- 
tial in the region between the reference plane and the charged surface. 

The nature and substance of the iterative schemes referred to in 
the above discussion are defined depending arbitrarily on the data that 
are available and the ingenuity of the programmer. A particular iter- 
ative scheme may have different convergence properties when applied to 
a number of different problems. A particular iterative scheme may have 
better convergence properties than other schemes when applied to a 
common problem. Furthermore, the selection of data defining a given 
problem may influence the convergence of a particular iterative scheme 
and the selection of an initial approximation may also influence the 
convergence toward a potential solution. Therefore, there are a broad 
range of topics that can be investigated and discussed relative to the 
synthesis of electric potentials by the iterative addition of ideal 
multipole fields. For the purposes of this report, the effects of 
boundary geometry, data selection and initial potential estimation will 
not be considered. Attention will be focused on the development of 
optimal iterative schemes based upon a well-defined model problem. 
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II. THE BOUNDARY VALUE PROBLEM 


The classical boimdary value problem in electrostatics may be 
stated as the problem of calculating values of a potential field exist- 
ing on a spatial region devoid of electric charge (i.e., where Laplace' 
equation is valid) when given the value of the potential on a simple 
boundary enclosing the region of interest. When the potential is 
specified as a continuous function on the boundary, the problem can be 
solved by a separation of variables technique and the expansion of the 
potential in terms of a complete set of orthogonal eigenfunctions 
(ref. 2). The geometric shape of the boundary Influences the choice 
of an appropriate coordinate system and the particular set of eigenfunc 
tions used for the expansion of the potential field. The solution is 
continuous inside the boundary and the error depends on the number of 
terms used in the eigenfunction expansion. 

In a practical situation, the value of the potential may be known 
only at a finite number of points on a boundary. The separation of 
variables and eigenfunction expansion method can be applied to obtain 
a continuous solution, if some approximation is used to interpolate a 
continuous functional value of the potential on the boundary from the 
finite number of known boundary values. The error of the solution de- 
pends on the approximation of the potential on the boundary as well as 
the number of terms in the expansion of the potential. The method of 
finite differences may also be used to approximate the potential inside 
a closed boimdary when only a finite number of boundary values are 
known (ref. 3). The solution is developed by using Taylor series ex- 
pansions to derive a finite number of linear potential equations in 
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terms of an equal number of unknowns (e.g., the potential values at all 
points on the grid enclosed by the boundary) and a finite number of 
knowns (e.g., the potential values at all points on the grid coincid- 
ing with the boundary) . A solution may then be obtained by using 
Gaussian elimination procedures to solve the system of linear equations. 
A finite difference solution is not continuous since it is defined only 
at nodal points on the grid. Therefore, interpolation errors will 
exist in addition to Taylor approximation errors, both of which depend 
on the spacing of the grid lines. Boundaries with irregular shapes will 
be problematic; since either a very fine grid will be required, leading 
to lengthy computer calculations, or some boundary approximations will 
be required to accommodate a larger grid spacing. 

A. Dipole Synthesis of Potential Fields 

Consider an iterative scheme for summing the potential fields of 
ideal electric dipoles to develop solutions to Laplace's equation on a 
domain enclosed by a boundary where the potential is known at a finite 
number of points (i.e. , Dirlchlet boundary condition). An initial esti- 
mate of the potential satisfying Laplace's equation is an average of 
the known boundary values plus any linear terms suggested by the bound- 
ary values (e.g., 4>^(x,y) = a + bx + cy) . This estimated potential can 
be improved by the iterative addition of an ideal dipole outside the 
boundary and nearest the location exhibiting the largest discrepancy 
between estimated and known boundary values. The potential field of an 
ideal dipole is itself a solution to Laplace's equation. Hence, the 
linear sum of the fields due to many such electric dipoles is also a 
solution to Laplace's equation. We will assume to have synthesized a 




particular solution to the boundary value problem when the sum of all 
the absolute discrepancies between estimated and known boundary poten- 
tial values is reduced below some tolerable limit. Such a solution is 
not unique, since the choice of dipole locations is arbitrary. If the 
above procedure for the reduction of discrepancy is pursued indefinite- 
ly, ultimately one dipole of optimal magnitude will be located outside 
the boundary near each one of the points where a boundary value is 
given. Thus, the iterative scheme expands the potential as a linear 
sum of the fields associated with dipoles located outside the 
boundary. If the dipoles are allowed Infinitely close to the boundary 
points, then each dipole magnitude will be chosen to reduce the dis- 
crepancy existing at one boundary location and have no effect on the 
discrepancies at the other boundary locations. Solutions developed in 
this manner are unacceptable since the field singularities introduced 
near the boundaiy produce only localized reduction of discrepancy. Con- 
sequently, the continuity of the potential solution on the boundary is 
jeopardized by placing dipoles too near the boundary itself. There- 
fore, we shall expect to synthesize solutions to the boundary value 
problem, if the dipoles are located well outside the bounded domain of 
interest. Convergence (e.g., discrepancy reduction) will be slower as 
dipoles are moved further from the boundary, but better approximations 
will result on the boundary between the pointwise specified values. 
Irregularly shaped boundaries can be easily accomodated by the iterative 
procedure for dipole synthesis. In addition, the estimated potential 
is defined in a continuous manner as the sum of dipole fields at all 
points on the bounded domain and no interpolation from a coordinate 
grid is required. 
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To define a specific class of boundary value problems, assume that 
the boundary S of region A in the two-dimensional r-9 plane can be rep- 
resented by N points distributed around the origin with no two points 

having the same 9 coordinate value. Assume also that the potential 

2 

\j)(r,d) satisfies Laplace's equation (i.e., V = 0) and is given at the 

same N points on the boundary 

iKR^.en) = ; n = 1,2,3---N (II-l) 

In order to allow calculations of the potential on the interior region 

A, as well as the boundary S; construct a 2N x N coordinate grid with 

its origin at the geometric center or centroid of region A and including 

2 

the N points of equation (II-l). Define the 2N points (r^, 9^) letting 
r = R for m = 1,2, 3,* • • N and r = for m = N+1, N+2,* • • 2N; 

where R is the radius of the smallest circle having its center at the 
origin of the coordinate grid which also encloses the region of interest. 

In order to begin an iterative procedure for approximating the 
potential at points on the interior, given the potential at N points on 
the boundary, define an initial estimated potential i})^ existing as an 
approximation to ip. 


4>„(r ,9„) = l/2[i(; + 1/2 r cos 0 (II-2) 

omn n n nnmn 


This choice is arbitrarily based on the maximum and minimum values of 
the N known boundary potentials; it satisfies Laplace's equation (l..e.. 


V 4> =0) and it guarantees that the difference between ip and (J> on the 

o o 

boundary is no larger than minus 

As a general procedure for improving the accuracy of the estimated 
potential on the boundary S, consider the iterative addition of dipole 
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potential fields ,0 ) due to the ideal dipoles located at 

d m n 

(r^,0^) outside of the boundary S enclosing region A. 


4>T^(r ,0 ) = (j) (r ,0 ) 
m n ^o m n 


E 

I <}>,(r ,e„) 

d=l ° 


(II-3) 


The following algorithm is used for calculating in two dimensions 

(i.e., the r-0 coordinate system) the potential fields associated with 

ideal electric dipoles at (r,,0.), 

d d 


A / a \ A..[r,- r cos(0 -0,)] 

(l).(r,0) = d'- d m n d' 

d m n —z z y 

2r.r cos (0-0.) 
d m dm n d 


(II-4) 


as derived in Appendix A. We will assume that the estimated potential 

is a good approximation to tp on the interior of region A, as well as on 

the boundary S, when a sufficient number of dipole potentials (}), have 

d 

been added to cf)^ such that matches the boundary potential at all 
points where the latter quantity is given. 


(f) (R ,0 ) = i|^(R ,0 ) ; n = 1,2,* --N 
£> n n n n 


(II-5) 


The arbitrary constant A, = qAr represents the charge magnitude q of an 

d 

ideal dipole having separation Ar as measured along a radius extending 

from the centroid of region A to the location (r, , 0,) of the dipole. 

d d 

For the purpose of calculating the error of the approximation of ip 

by 4> , define the discrepancy function D as the difference between ijj 

and ^ at each of the N boundary locations where ip is known. 

E 


D 

n 


= 'I'CR ,0„) - 4>„(R ,9„) 

n n E n n 


(II-6) 


The average discrepancy and average absolute discrepancy, as defined 
by equations (II-7) and (II-8), -can be used as measures of the accuracy 
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of the estimated potential. 


D 


1 

N 


N 

I 


n=l 


D 

n 


(H-7) 


W - i ! IDJ 

n=l 


(n-8) 


We are interested in developing iterative procedures for reducing D 

and |d 1 . A uniform shift of the estimated potential can be used to 

reduce D to zero. Suppose the average discrepancy D exists for 

estimated potential (() . If a new estimated potential (()' is formed by 

adding the average discrepancy D to (j> , then we wish to show that the 

£ 


average discrepancy D for ({)' is equal to zero. 

£ 

, , N T N 

“ ’ i z J 

n=l n=l 


■ » \ 
n=l 


N 


(II-9) 


Then by rearranging terms and recalling the definition of D expressed 
by equation (II-7) , we derive the desired result. 

“ I ", i 15-5-5 = 0 (ii-io) 

n=l n=l 

Thus, the estimated potential can be shifted so that its average 
discrepancy with known boundary values is reduced to zero. 

The reduction of the average absolute discrepancy is not as 
easily accomplished. Consider the rotation of the estimated potential, 
with respect to the coordinate grid, as one operation for reducing 
I'oy on the boundary. Define N rotated potentials 


R = 1,2,* '‘N such that 
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^ = 1.2,--«N-R (II-ll) 

'f'R^’^m’V = ‘<>^’^m’VR-N> ^ = N-R+1. • • -N (H-12) 

for all m = 1,2, •••2N. Then one might choose, as the best estimated 
potential, the rotated potential yielding the lowest average absolute 
discrepancy and zero average discrepancy. 

The average absolute discrepancy might also be reduced by scaling 
the estimated potential after the rotating and shifting operations. 

The following algorithm is proposed for scaling the estimated potential 
so as to reduce the average absolute discrepancy whenever the average 
discrepancy is zero. 

N 

I 

I . 

n=l 


♦ ' i *«n-V 

n=l 


(11-14) 


n— 1 


(11-15) 


After the rotation, shifting and scaling operations, the addition of 
ideal dipole potentials can be used to further reduce |d| and improve 
the approximation of by (}). The definition of a specific iterative 
scheme depends on the general criteria that can be established for 
determining optimal locations and magnitudes for each successive 
electric dipole that is added near the boundary enclosing the re- 


gion of interest. 



13 


B. Iterative Schemes for the Dipole Synthesis of Solutions to the 
Boundary Value Problem 

VM0DGR.F4 is the name of the main Fortran routine used for solving 
boundary value problems. Initially, subroutine DATA is called to set 
up the problem in the X-Y coordinate system corresponding to the user's 
particular reference frame. The subroutine asks the user for up to 
forty points representing the boundary and for the potential value at 
each of the points which must be specified in a clockwise sequence 
around any plane boundary. The main routine generates an r-0 coor- 
dinate grid with origin at the centroid of the region enclosed by the 
boundary . 

The determination of a lower bound on the distance between the 
boundary and the dipole locations is the primary objective of the 
study. It is assumed that all dipoles will be best aligned if their 
axes, extending in the direction of charge separation, pass through the 
centroid of the bounded domain. Subroutine ADPOLE locates a dipole at 
the angular coordinate 9^ where the discrepancy D is greatest between 
the estimated potential (|) and the boundary potential ip. Then a param- 
eter G is defined in subroutine ADPOLE so that all dipoles are located 
a fixed radial distance from the centroid. This distance 

r^ = R(l+G/N) (11-16) 

is calculated in terms of the number N of known boundary potentials 
and the largest radius R needed to construct a circle centered at the 
centroid and enclosing all of the region of interest. Dipoles will be 
allowed closer to the boundary as the value of the G parameter, speci- 
fied by the user interactively with subroutine DATA, is allowed to 
decrease toward a limiting value of zero. Once a dipole location has 
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been determined, subroutine ADPOLE choses a dipole magnitude A, that 

d 

reduces the discrepancy at the boundary location nearest the dipole 

to approximately zero. Alternatively, subroutine TADPOLE may be called 

by the main routine to select dipole locations and magnitudes. In this 

subroutine, dipoles are also located at the angular coordinate 9, 

d 

corresponding to the greatest boundary discrepancy and the dipole 
axis is aligned with the centroid of the region of interest. However, 
for the TADPOLE scheme, the G parameter does not fix the radial coor- 
dinate for dipole location but specifies a minimum value that can be 
allowed. Whenever possible, a dipole magnitude A^ and dipole radius 
r^ are chosen simultaneously so that 

' -“l 

and 

V “ ■ IdJ ( 11 - 18 ) 

where the coordinates (R ,9 ) represent the point on the boundary nearest 

c c 

the point where the discrepancy has a sign opposite from the sign 

of the largest discrepancy D^. Whenever, the criterion expressed in 
equations (11-17) and (11-18) cannot be achieved without placing the 
dipole inside the minimum radius, as expressed by equation (11-16), r^ is 
chosen equal to this minimum radius on the basis of the G parameter 
specified in subroutine DATA. Thus, the T parameter, also entered 
interactively with subroutine DATA, balances discrepancy reduction at 
the boundary point where the discrepancy is maximum with the accompaning 
discrepancy increase at a nearby boundary location. 

The main routine VM0DGR.F4 first defines an initial estimated 
potential according to equation (II-2) and then calls subroutine 
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ROTATE to achieve discrepancy reduction by the rotation of the initial 
estimated potential relative to the coordinate grid. An iterative 
scheme, calling subroutines SHIFT, SCALE and either ADPOLE or TADPOLE; 
is then applied to accomplish the reduction of boundary discrepancies. 
Subroutines SHIFT and SCALE, listed in Appendix B, accomplish discrep- 
ancy reduction according to equations (II-IO) and (11-13) respectively. 
The other Fortran routines described in this section are also listed 
in Appendix B. 

The particular problem of a square boundary with a potential of 
10 volts specified on three sides and 0 volts specified on the fourth 
side has been chosen to demonstrate the method of dipole synthesis and 
to Investigate various criteria for choosing dipole locations in an 
iterative scheme. Figure 2 shows how twelve boundary points were 
selected at uniform angular intervals and how one interior point at 
(r ® .471, 9 = 90®) was selected for comparing synthesized potentials 
with a precise series expansion which yields a value of 9.59 volts 
just inside the unit square boundary. Tables 1 and 2 show the poten- 
tial values synthesized at the test point with the main routine VMODGR 
using subroutines ADPOLE and TADPOLE, respectively. The fields of 
ideal dipoles were added at each boundary point and at the test point 
in an iterative process, one dipole per iteration, until the sum of the 
absolute discrepancies was reduced below 10 ^ volts. Note, in Table 1, 
that the estimated potential converges to its precise value at the test 
point as the dipoles are moved further from the boundary. However, the 
number of iterative computer calculations required is large. When di- 
poles are allowed closer to the boundary (i.e., smaller G parameter), 
convergence is much faster but errors in the estimated potential at the 
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Table 1 


Potential Synthesis With Constant Dipole Radius 


G 

Iterations 

Test Potential 
(volts) 

1 

65 

10.42 

2 

73 

10.08 

4 

109 

9.83 

2tt 

157 

9.67 

10 

354 

9.61 

15 

791 

9.59 

20 

1718 

9.59 


Table 2 

Potential Synthesis With Variable Dipole Radius, G = 2ir 
T Iterations Test Potential 




(volts) 

0.01 

157 

9.67 

0.30 

157 

9.67 

0.50 

161 

9.64 

0.70 

169 

9.59 


0.90 


214 


9.59 
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test point are substantial. Note, in Table 2, that the G parameter 
is set equal to 2iT thereby establishing the minimum radial coordinate 
available for the location of point dipoles. The convergence of the 
estimated potential to the specified boundary values is fastest when 
the T parameter forces many dipoles close to the boundary (i.e., smaller 
values of T) . The occasional placement of dipoles at greater distances 
from the boundary, using both the G and T parameters, provides a flex- 
ibility that allows compromise between the speed of the convergence and 
the accuracy of the estimated potential at interior locations. Thus, 
convergence with the required accuracy at the test point can be achieved 
with 791 iterations using the G factor alone and with 169 iterations 
using the G and T factors together. 

Once convergence is achieved, VMODGR calls subroutine VGRAPH to 
plot equipotential lines as determined on the interior coordinate grid. 
Figures 3 and 4, corresponding to G parameter values of 2 and 2tt re- 
spectively, show solutions to the boundary value problem for 20 points 
on the square boundary spaced at equal angular intervals. Figure 5 
shows the potential solution for a G parameter of I'n and a T parameter 
value of 0.1. A comparison of Figures 4 and 5 reveals no major dif- 
ferences between the potential solutions generated by the two iterative 
schemes. However, the solution shown in Figure 5 was generated by 
using subroutine TADFOLE and required less than half the computer time 
needed when subroutine ADPOLE was used. Figure 6 shows the improvement 
in the solution potential resulting from the specification of 36 bound- 
ary values as compared to only 20 boundary values. It should be noted 
that the grid structure used Inside the boundary results in linear 
interpolation errors in the drawing of equipotential lines. These 
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errors are not a consequence of the procedure for potential field 
synthesis but result when an attempt is made to display a potential 
field based upon its value at a finite number of points. The potential 
display could be improved by a better interpolation procedure or by 
calculating the potential at more points. However, the accuracy of 
the potential estimate synthesized at a particular point of interest 
does not depend on the number of interior points considered. Therefore, 
potential estimates could be synthesized at as many, or as few, in- 
terior points as desired with accuracy depending on the number of 
boundary values specified. 



24 


III. THE DEFINITION OF A MODEL PROBLEM 


The electric potential resulting from a distribution of surface 
charge can be shown to satisfy Laplace's equation everywhere except 
on the surface itself. Since the potential fields associated with 
ideal multipoles are also solutions of Laplace's equation, we assume 
that a linear sum of weighted multipole potentials can be used to 
construct a particular solution that approximates the potential field 
created by the surface charge distribution. This particular solution 
will have singularities at the spatial locations of the multipoles. 
Therefore, the spatial domain near a charged surface can be partitioned 
into two subspaces; one subspace defined as the region where a valid 
potential approximation is desired and another subspace reserved for 
the location of ideal multipoles. In the boundary value problem con- 
sidered previously, the closed boundary separated the interior subspace 
(where solutions were desired) from the exterior subspace (where dipoles 
were located) . In a general situation, the charged surface may only 
partially define the subspace where multipole location is allowed. In 
any event, an iterative scheme for summing multipole potentials to 
approximate the potential near a charged surface will depend on the 
geometry of the particular problem under consideration. 

Furthermore, in contrast to the boundary value problem, the sur- 
face charge and surface potential may be unknown in the general problem. 
Trajectory data for electrons moving in an unknown potential field 
created by an unknown charge distribution can be used as constraints 
on the class of functions that approximate the potential field (ref. 4). 
Electron deflection measurements can be defined, in the absence of 
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boundary values, as a basis for calculating the discrepancy between 
the potential in the region probed by the electron trajectories and 
the multipole potential approximation generated by an iterative pro- 
cedure. In such a case, the electron trajectories define a subspace 
where the potential approximation is desired and where ideal multipoles 
cannot be located. Therefore, a particular iterative scheme for deter- 
mining the optimal locations of dipoles, quadrupoles or multipoles will 
be influenced by the electron trajectory data provided as part of a 
general problem. 

To define a particular problem that can be used to compare the 
convergence properties of various iterative procedures, consider the 
potential field created by two uniformly biased fins in the vicinity 
of a zero potential ground plane. So that calculations may be per- 
formed in a two-dimensional space, assume that both the ground plane 
and the biased fins extend infinitely in the positive and negative z 

coordinate directions. Note, in Figure 7, that the fins begin at x-y 

DX 3DX 

coordinate values (-^,DY) and (— ^,DY) respectively and extend in- 
finitely in the positive y dimension. Suppose the x-y coordinate plane 

can be mapped onto the x'-y' coordinate plane by shifting the origin by 

DX DX DY 

amounts DY and — , by scaling the axes by factors of — and — , and by 

rotation of the axes through 90®. 


x' = ^ (DY-y) and 7* = ^ (x-|^) 


(III-l) 


The conformal mapping of the complex z plane (z = x' + iy*) onto the 
complex w plane (w = u + iv) is 


w = sin ^(e^) 


(III-2) 
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which collapses the infinite fins onto a finite interval (i.e. , [0'',0’]) 
along the u-coordinate axis (ref. 5). The value of this transformation 
lies in the fact that the potential field created by the biased fins 
in the x-y plane transforms to a potential with a uniform gradient in 
the u-v plane. Thus, 


(j)(z)-C = v(z) = Im[w] = Im[sin ^(e^)] 


(III-3) 


“1 z 

can be used with an appropriate expression for sin (e ) to construct 


an algorithm 


, 1/2 


<{)(x',y')-C = Im{-i*ln[ie *e ^ + (1-e *e ^ ) ]} (III-4) 


that can be used to calculate the biased fin potential at selected 
points in the x-y plane. The constant C is required so that the 
potential can be shifted enough to simulate the zero potential line 
X* = DY (i.e., y = 0) . In order to simplify the algebraic evaluation 
of the natural logarithm in equation (III-4) , convert the complex 
argument of the square root to polar form. 

1/2 


r = (1+e -2e^ »cos 2y') 


(IH-5) 


1 2x' . , ’ 

- _ . ^ -Ir-e 'Sin , 

= TT + tan [ ^ — J . 

T 2x' , , 

1-e 'COS ^y 


(III-6) 


These expressions lead to the form 


(j)(x’ ,y' , r, 9) = Im[-i *ln(ie^'e^^ + /re^^^^)] + C. (III-7) 


To further simplify the calculations, let 


R = [r+e^^ -2»^e^'sin(y'-0/2) 


(III-8) 



/ 
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and 


, x' t . j 0/2 

Q = IT + tan'^^ e cosy + /F sin ^ 

*/r cos 6/2 -e^ sin y' 


(HI-9) 


to obtain 


4)(R)-C = Im[-i*ln(Re^)] = Im[-i*ln(R)±J^] = -ln(R) . 


(III-IO) 


Therefore, to find the potential anywhere in the x-y plane; we first 
calculate values of x' and y', then find the corresponding values of 
r and 9, and finally determine a value of R. The potential at point 
(x,y) is then equal to the natural logarithm of R plus an arbitrary 
constant . 

This algorithm for calculating the potential field in the x-y 
coordinate space of the biased fins is implemented by main routine 
TEST.F4 (ref. Appendix C) . When the additive constant is properly 
chosen, the potential falls to zero along the x axis within a tol- 
erance of 1% of the fin potential if the spacing of the fins is not 
greater than the distance from the fins to the x axis by more than a 
factor of 1.8. Thus, the gound reference plane could be moved to 
within a distance DY = DX/3.6, when the spacing' of the fins is DX/2, 
without altering the potential field of the biased fins by significant 
proportions. TEST.F4 calculates potential values on an 80 x 27 
coordinate grid corresponding to a region DX by 1.4DY in the x-y plane. 
The I and J indices are used to identify points in a region with dimen- 
sions DX by DY according to 


(III-ll) 


and 


X(I) = (I-2)DX/80 
Y(I) = (J-2)DY/20. 


(III-12) 
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These values are written into a storage file POTL.DAT which may be 
accessed by other routines that require the test potential. Table 3 
lists values of the biased fin potential at 110 points within the 
spatial regtion of interest. TEST.F4 also includes a graphic sub- 
routine, VPLOT, that utilizes a linear interpolation technique to 
plot equipotential lines in the space between the biased fins and the 
ground reference plane. Figure 8 shows the graphical plot generated 
by TEST.F4 with input parameters, DX = 3.6 and DY = 1.0. 

Also consider six electron trajectories, emanating at six loca- 
tions along the x axis, that probe the region between the fins before 
falling back to the reference plane. Assume that the initial coordinate 
values and the initial coordinate velocities are known as well as the 
final coordinate values of the six trajectories that serve as con- 
straints placed upon approximations of the biased fin potential. In 
particular, let five trajectories emanate from locations (0.900,0), 
(1.305,0) , (1.710,0), (2.105,0) and (2.520,0) along the x axis and 
with an initial direction parallel to the y axis and with initial 
kinetic energy equal to 85% of the fin bias. Also let a scaling tra- 
jectory emanate from the location (0.855,0) on the x axis with an 
initial velocity in the x direction equal to 1.1 times the velocity 
in the y direction. The main routine DATA.F4 is used to determine the 
final coordinates of the six trajectories defined by the initial electron 
parameters and the test potential. DATA.F4 calls subroutine DVOGEL in 
an iterative manner to calculate successive electron steps through a 
test potential specified by reading file POTL.DAT. DVOGEL solves the 
electron's second order equations of motion according to a method pro- 
posed by the mathematician DeVogelaire (ref. 6). Electric fields 
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Table 3 


Values of the Test Potential at 110 Selected Points 
y-coordinate index (J) 


21 

2.02 

1.82 

1.67 

1.57 

1.52 

1.50 

1.52 

1.57 

1.67 

1.82 

2.02 

19 

1.69 

1.61 

1.50 

1.43 

1.38 

1.37 

1.38 

1.43 

1.50 

1.61 

1.69 

17 

1.44 

1.40 

1.33 

1.27 

1.24 

1.23 

1.24 

1.27 

1.33 

1.40 

1.44 

15 

1.23 

1.20 

1.16 

1.11 

1.09 

1.08 

1.09 

1.11 

1.16 

1.20 

1.23 

13 

1.02 

1.01 

0.98 

0.95 

0.93 

0.92 

0.93 

0.95 

0.98 

1.01 

1.02 

11 

0.83 

0.82 

0.80 

0.78 

0.76 

0.76 

0.76 

0.78 

0.80 

0.82 

0.83 

9 

0.65 

0.64 

0.63 

0.61 

0.60 

0.59 

0.60 

0.61 

0.63 

0.64 

0.65 

7 

0.46 

0.46 

0.45 

0.44 

0.43 

0.43 

0.43 

0.44 

0.45 

0.46 

0.46 

5 

0.28 

0.28 

0.27 

0.27 

0.26 

0.26 

0.26 

0.27 

0.27 

0.28 

0.28 

3 

0.10 

0.10 

0.10 

0.09 

0.09 

0.09 

0.09 

0.09 

0.10 

0.10 

0.10 





x-coordinate 

index 

( I ) 
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Figure 8. Equipotential Lines for the Test- Biased Fin 
Potential on The Domain 22<I<62, 2< J <22 
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required by subroutine DVOGEL at points between the grid lines are ob- 
tained by a call to subroutine NEWTON which implements a divided 
difference interpolation (ref. 7) of the electric fields based upon the 
grid potential values provided by POTL.DAT. The electric fields at the 
16 grid points (i.e., a 4 x 4 minigrid) closest to the point where an 
interpolated electric field is required, are used to determine a spatial 
step size (not in excess of 0.01 times a typical diagonal measurement of 
the region of interest) used for the piecewise development of electron 
traj ectories. 

The object of routine DATA.F4 is the determination of the final 
coordinate values (0.901,0), (1.848,0), (1.824,0), (1.709,0) (2.086,0) 
and (2.691,0) for the six constraining trajectories. These final 
coordinate values, in addition to the initial coordinate values and 
velocities are written into file TRAJ.DAT which can be referenced by 
routines that approximate the potential field of the biased fins. 

DATA.F4 includes a graphics routine for plotting the constraining 
electron trajectories and calls subroutine VPLOT to superimpose 
equipotential lines representing the test potential as shown in 
Figure 9. 



Figure 9. Equipotential Lines for the Test-Biased Fin Potential 
and Constraining Electron Trajectories on the Domain 2< I <82, 2< J <22 


IV. DIPOLE SYNTHESIS OF ELECTROSTATIC POTENTIALS 


Preliminary investigation was begun with the intention of de- 
veloping methods for dipole synthesis of three-dimensional electrostatic 
potentials assuming a uniform potential in one coordinate direction. 

The problem can be considered two dimensional if a line dipole is con- 
structed by letting two infinite lines of opposite charge approach each 
other in the limit of superposition. The result is a line dipole 
potential which is uniform in the coordinate direction parallel to the 
line dipole. Furthermore, in a plane perpendicular to the line dipole; 
a difference results from integration along the line charges pro- 
ducing a 1/r variation of the potential with distance from the line 

2 

dipole, rather than the 1/r dependence associated with a point dipole. 

The line dipole described above is to be used as the basic build- 
ing block for synthesizing unknown potentials in a two-dimensional 
space that can be probed by electron beams. This space is bounded 
above by the charged surface that produces the potential of interest 
and bounded below by the zero potential reference plane of the electron 
beam apparatus. Since the potential must be zero on this boundary, 
line dipoles can be placed in pairs, one on either side of the boundary 
so as to be electrostatic images of each other. Figure 1 shows the 
space of interest, the zero-potential boundary and a line dipole pair. 
Notice that the pair of line dipoles is located outside the space of 
interest so that Laplace's equation will hold on the interior. 

The model problem under consideration requires the determination 
of the magnitude, location and number of line dipole pairs that best 
approximate the potential existing in the space of interest resulting 



from an unknown surface charge distribution. The potential in the 
region of interest is described by electron trajectory data including 
initial electron coordinates and velocities and final electron coor- 
dinates. The basic method begins with a uniform field approximation 
of the potential which can be modified by the iterative addition 
of line dipole potentials until a potential is created that reproduces 
the electron trajectory data. 

To formulate an iterative scheme for the dipole synthesis of 
solutions to the model problem, suppose that the potential ij; exists 
as a result of the biasing of two fins relative to a zero potential 
reference plane. The last of the six trajectories listed in the state- 
ment of the model problem (i.e., the scaling trajectory) can be used 
for a uniform field approximation of tp in order to develop an initial 
estimated potential (J)^. The equations of motion describing the scaling 
trajectory are 



3(f) 

= -5—^ = 0 and 

m 0X 




m 9y m 


A . 

o 


(IV-1) 


A linear system of two equations 


XT(6) = X0(6) + T'VX0(6) (IV-2) 

and 

YT(6) = Y0(6) + T-VY0(6) + J ^ A • T^ (IV-3) 

z m o 

in two unknowns (i.e., the constants A^ and'T) is developed by integrating 
equations (IV-1) twice and substituting the data for the scaling tra- 
jectory. When YT(6) = Y0(6) = 0, the result is 
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_ ^ VX0(6)-VY0(6) 
o ~ q [XT(6)-X0(6)]* 


(IV-4) 


Thus, an initial estimated potential, d) =A y, can be calculated at all 

o o 

points in the region of interest based upon the scaling trajectory data 
and the electron charge to mass ratio The estimated potential can 
be synthesized from (j) and a number of dipole fields (see Appendix A) , 


(j>^(x,y) = A^[ 


y-yd 

(x-x^)^ + (y-y^)^ 


y^d 

(x-x^)^ + (y+y^)^ 


]. 


(iv-5) 


to improve the estimation of the biased fin potential, 

c{)(x,y) = (l)^(x,y) + Z (f>^(x,y), (IV-6) 


according to some iterative scheme that specifies optimal dipole param- 
eters A,, X., y, for meeting the constraints imposed by the electron 
d d d 

trajectory data (i.e., trajectories 1-5). For each iteration of the 
scheme, estimated electron trajectories are developed from the initial 
trajectory data (i.e., X0(K), Y0(K); K=l,2,...6) and the equations of 
motion (IV-1) substituting, for (j)^, the best estimate of the potential 
field (})(x,y) calculated from equation (IV-6). The coordinates XM(K) 
and YM(K) at the highest points on the estimated trajectories are re- 
corded as well as the exit coordinates XE(K) and YE(K); for K=l,2,--*6. 
A discrepancy is defined 


D(K) = 


rXT(K)-XE(K)1 

[XT(K)-XO(K)] 


K=l,2,---5 


(IV-7) 


based upon estimated trajectory data XE(K) and actual trajectory data 
XT(K) and XO(K). In order to reduce the discrepancies, dipoles should 
be located between trajectories having discrepancies opposite in sign. 
Therefore, the trajectory with the largest absolute discrepancy 
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(e.g., K=K^) is identified and the nearest trajectory having a dis- 
crepancy of opposite sign (e.g. , K=K.) is identified. The dipole x, 

2 . d 

parameter, 

j [XM(K^) + XMCK^) 1 , (IV-8) 

is chosen midway between these two trajectories. The dipole y, param- 

d 

eter is fixed, y^=1.4-DY, at a reasonable distance from the region of 
interest (e.g. , as close to the region as possible without introducing 
local singularities, as discussed for the boundary value problem). 

Then a dipole magnitude can be chosen which reduces the discrepancy 
of the trajectory to approximately zero. 

Toward this end, consider the integration of the equations of 
motion (IV-1) , with 4) replacing (p^, for electrons following the 
estimated trajectory through the estimated potential. If the time T 
represents the time required for electrons to traverse this trajectory, 
then the results of the integration are 


XE(K^) 


X0(KJ + T-VX0(KJ + ^ 

1 1 / m 


'0 ■'( 




and 


YE(K^) 


Y0(K,) + T-VY0(K,) 

1 1 Z m 


rt’ 


9y 


dt dt’ 


(IV-9) 


(IV-IO) 


The trajectory consists of a number of points whose coordinates have 
been determined by a stepwise solution of the equations of motion. The 
electric fields, and -|^, at these points were calculated along with 
the time steps for electron motion between the points in the development 


of the trajectory. Therefore, the integrals in equations (IV-9) and 


(IV-10) can be replaced by a pointwise summation of the electric fields 
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over the N points of the trajectory using trapezoidal integration to 
develop 


XI (Kj^) = 


rT 

■0 


ft' 

■'o 


3x 


dt dt' 



(S 


n-1 


+ S )At 
n n 


where 


(IV-11) 


and 



(IV-12) 




fT 

0 


Jq 3y 


dt dt' 



(R 


n-1 


+ R )At 
n n 


where 


(IV-13) 


R 

n 


1 

2 


<^yi_i 


A 


i=l 


+ 

% 


)]At . 
i ^ 


(IV- 14) 


Substituting the approximation of equation (IV-13) into equation 
(IV-10), we have for YE(Kj^) = YO(Kj^) = 0, 


T = 


2 m VYO(K^) • 


(IV-15) 


Substituting the approximation of equation (IV-11) into equation (IV-9) , 
we have 

XE(K^) = XO(K^) + T-VXO(K^) +j^XI(Kj^). (IV- 16) 


Now, we wish to add a dipole potential 4>^ to the estimated potential 
<j), altering the trajectory only slightly and shifting the exit 
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coordinate from XE(K^^) to the data coordinate XT(K^). The equations 
of motion for this new estimated potential are 


= £ /M + 

m '■3x dx 


and 


^ = 5. /li 
at^ m ^9y 



(IV-17) 


Integration of the first of these equations (IV-17), yields 


XT(Kj^) = XO(K^) + T-VXO(IC|^) 



fT 

0 





3(j) 

^)dt dt- . 


(IV-18) 


When trajectory is shifted only slightly, the integral in equation 
(IV-18) can be replaced by the summation of electric fields over the 
same N points used to generate equation (IV-16) . Hence, 


XT(Kj^) = XO(K^) + T-VX0(Kj_) + [XI(K^) + A^-VI(Kj_)], (IV-19) 

where 

1 ^ 

VI(K.) = j E (V + V )At (IV-20) 

1 / 1 n-1 n n 

n=l 

and where 


V 

n 




•) 

i-1 


+ ( 


d<P 


(IV- 21) 


Subtracting equation (IV-16) from equation (IV-20) , we may solve for 
the dipole magnitude 


^ XT(K^) - XE(K^) 
""d q ^ VI(K^) ^ 


(IV-22) 


After the addition of a dipole with magnitude A, at coordinates 

a 

(Xy»y^), one iteration of the scheme for discrepancy reduction is 
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completed by scaling the dipole modified estimated potential, 

<{>' (x,y) = SF*4)(x,y) . 

The scale factor, 

^ XE(6) - X0(6) 

XT(6) - X0(6)’ 

is calculated using the exit coordinate XE(6) obtained for the stepwise 
development of the scaling trajectory in the dipole modified estimated 
potential. 

Appendix D lists the main Fortran routine PTENT.F4 which imple- 
ments the iterative routine described by equations (IV-1) through 
(IV-24) . The user has the option of entering electron trajectory data 
by reading from file TRAJ.DAT or by answering pertinent questions in an 
interactive user mode. Table 4 shows the reduction of discrepancy and 
and of error estimates obtained for 15 iterations of the routine. Note 
that these quantities do not decrease monotonically and limits are 
reached after about 8 iterations of the routine. Figure 10 shows 
equipotential lines drawn with a linear interpolation subroutine of 
PTENT similar to the VPLOT subroutine described for the boundary value 
problem. The solution to the model problem depicted in Figure 10 can 
be compared to the test potential in Figure 8. 

Appendix E contains a listing of the main Fortran routine RTENT.F4 
which is also an iterative routine for the dipole synthesis of solutions 
to the model problem. RTENT.F4 differs from PTENT. F4 only in the 

methods used for determining the dipole x, ,y, and A factors. For this 

d d d 

new scheme, the x, parameter can be calculated by equation (IV-8) when 
d 

is identified as the trajectory having the largest positive dis- 
crepancy and K 2 is identified as the trajectory having the largest 


(IV- 23) 


(IV-24) 
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Table 4 

Error Estimates for Dipole Synthesis with PTENT.F4 


Iteration 

Discrepancy ,DD 

Average 

Error 

Average 

Absolute 

Error 

Maximum 

Absolute 

Error 


5 

E lXT(K.) - 

XE(K^ 1 

(%) 

(%) 

(%) 

0 

i=l ^ 

1.51 


4.6 

5.8 

20.0 

1 

1.29 


3.5 

5.0 

21.6 

2 

1.28 


1.9 

4.6 

20.0 

3 

1.33 


1.7 

3.8 

16.8 

4 

0.92 


3.7 

4.4 

18.3 

5 

0.92 


3.0 

4.2 

18.1 

6 

0.70 


2.3 

2.6 

14.4 

7 

0.45 


2.9 

3.7 

17.9 

8 

0.42 


1.7 

3.5 

16.8 

9 

0.45 


3.0 

3.7 

17.9 

10 

0.43 


1.8 

3.4 

16.8 

11 

0.45 


3.0 

3.7 

17.9 

12 

0.43 


1.9 

3.4 

16.9 

13 

0.45 


3.1 

3.7 

17.9 

14 

0.43 


2.0 

3.4 

16.9 

15 

0.45 


3.1 

3.8 

18.0 




Figure 11. Solution to the Model Problem Synthesized with RTENT.F4 
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negative discrepancy. The parameter is arrived at by a trial and 
error procedure designed to reduce the discrepancies of both the 
and trajectories. This is accomplished by finding a minimum value 
of 


3<f>d 

DCKf) ^[XM(Kj^),YM(Kj^)l 
D(K ) “ 

-^[XM(K2),YM(K2)] 


(IV- 25) 


as y^ is allowed to run through a range of values. When the values 

of X, and y , have been chosen, a value of A, is chosen to reduce the 
d d d 

discrepancy of the trajectory (or the trajectory if its dis- 
crepancy is larger in absolute magnitude than the trajectory) to 
approximately zero as described by equations (IV-9) through (IV-22) . 
Table 5 lists the error estimates for 15 iterations of the RTENT.F4 
routine as compared to Table 4. Note that error estimates are larger 
for RTENT.F4 than for PTENT.F4 and that the reduction of discrepancy 
is slower. However, RTENT.F4 consistently underestimates the biased 
fin potential while PTENT.F4 consistently yields an overestimation. 
Therefore, upper and lower bounds can be established for the solution 
of the model problem. Figure 11 shows the estimated potential gen- 
erated by RTENT.F4 which can be compared to Figures 10 and 8 to 
appreciate the accuracy of these iterative dipole routines. 

Tables 6, 7 and 8 show the potential error estimates at 110 
selected locations in the region of interest for 5, 10 and 15 iter- 
ations of routine RTENT.F4; respectively. First note that the largest 
errors occur at the edges of the region of interest in spatial regions 
unprobed by the electron trajectories. Note also that the errors at 
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for n- 

Oxpoie 

^>^“thesls vith 

RTBNT.P4 
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nearly all of the selected locations decrease monotonically from the 
fifth to the tenth to the fifteenth iteration. 
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V. QUADRUPOLE SYNTHESIS OF ELECTROSTATIC POTENTIALS 

The two iterative routines for dipole field synthesis of electric 
potentials yield good approximations of the biased fin potential. However, 
the addition of dipole fields individually chosen to reduce trajectory dis- 
crepancies may result in slow convergence and lengthy computer runs (e.g., 
the 15 iterations appearing in Table 4 required 30' 38" of CPU time). 

Efforts to speed the convergence of the synthesis process by adding dipole 
arrays, rather than single dipoles, have been unsuccessful. Note, in Fig- 
ure 12, that the x-component of the electric field associated with the line 
dipole is not symmetric about the axis of the line dipole. The five 
curves of Figure 12 represent the electric field in the x-y reference plane 
of the biased fins (see Figure 7 ) along a line midway between the reference 
plane (y=0) and the biased fins (y=DY), as the result of a line dipole and 
its image located at five perpendicular distances from the reference plane. 
Note also that the dipole has equal and opposite effects at two spatially 
distinct locations in the region of interest. These characteristics of 
the dipole field are problematic when considering the placement of several 
line dipoles simultaneously in an effort to reduce the discrepancies associ- 
ated with several electron trajectories. Consequently, the optimal place- 
ments of dipole array elements are not easily deduced since it is difficult 
to associate discrepancy reduction for one trajectory with a particular 
element of the dipole array. 

A. Quadrupole Synthesis of Solutions to the Model Problem 

Although the difficulties associated with synthesizing potential fields 
with dipole arrays may not be insurmountable, we have chosen to investigate 




Figure 12. Electric Fields for Five Dipole and Quadrupole 
Locations along a Line inside the Region of Interest 
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the quadrupole synthesis of electric potentials. In analogy to the 
ideal line dipole, consider the ideal line quadrupole consisting of four 
lines of alternating charge extending infinitely in a direction perpen- 
dicular to the two dimensional region of interest. Figure 12 shows 
these four lines of charge which are allowed to approach each other in 
the limit of superposition to derive the potential field of the ideal 
line quadrupole. Since, we are interested in the synthesis of poten- 
tials which fall to zero along a reference plane, an image line 
quadrupole is simultaneously located an equal distance from the refer- 
ence plane but on the side opposite the location of the line quadrupole 
itself. The potential field, as derived in Appendix A, due to a line 
dipole at (x ,y ) and its image at (x ,-y ) is 

q q q q 


(x,y)=A (x-x ) 

q q q 


(y-Hy_^ 




_[(x-Xq)^ + (y+y^)^]^ + [(X-X^)^ + (y-y^)^]^ J . (V-1) 


Notice that this potential vanishes along the reference plane (y=0) 
regardless of the values of the quadrupole parameters A^, x^, and y^. 

In analogy with dipole synthesis, the line quadrupole potential can be 
used for the synthesis of unknown potentials existing in a two dimen- 
sional space of the biased fins (i.e. , along a line, y=constant, midway 
between the fins and the reference plane) as the result of an ideal 
line quadrupole and its image located at the same five perpendicular 
distances that were considered for the ideal dipole field of Figure 12. 
In comparing these two diagrams, notice that the quadrupole field is 
symmetric about its axis. Furthermore, the quadrupole field is largest 
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on its axis, rather than at locations on either side of the axis as was 
found for the dipole field. These quadrupole field characteristics 
facilitate the definition of iterative schemes for the synthesis of un- 
known potentials. In particular, a number of ideal quadrupoles can be 
located simultaneously with that quadrupole which is placed directly 
over each trajectory being primarily responsible for the reduction of 
discrepancy associated with that trajectory. However, before consider- 
ing the topic of potential synthesis using the fields of quadrupole 
arrays, the synthesis of electric potentials by the Iterative addition 
of single quadrupole fields will be considered. This will provide a 
basis for comparing dipole and quadrupole methods for the synthesis of 
solutions to Laplace's equation subject to the constraints imposed by 
electron trajectory data. 

To formulate an iterative scheme for quadrupole synthesis of 
solutions to the model problem, suppose that the potential \p exists 
as the result of the biasing of the fins relative to a zero potential 
reference plane. An initial estimated potential (j)^ can be defined as 
the uniform-field approximation of ^ which satisfies the data for the 
last of the six trajectories listed in the statement of the model prob- 
lem (i.e., the scaling trajectory). In analogy with equations IV-1) to 
(IV-4) , an initial estimated potential = A^y, can be calculated at 
any point in the region of interest based upon the data for the scaling 
trajectory and the electron charge to mass ratio. After determining the 
uniform field approximation of the estimated potential, the estimated 
potential can be modified according to an iterative scheme that adds 
quadrupole fields 4>^(x,y) to cj)^(x,y) in order to improve the approxi- 
mation of the biased fin potential ij;(x,y). 
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Q 

(j)(x,y) = 4)g(x,y) <|)q(x,y) (V-2) 

For each iteration of the scheme, estimated electron trajectories in 
the estimated potential field are developed from a knowledge of initial 
electron coordinates and velocities for the first five electron tra- 
jectories listed as data for the model problem. The coordinates (XM(K) , 
YM(K)) at the highest points on the estimated electron trajectories are 
recorded as well as the exit coordinates (XE(K), YE(K)) for K = 1, 2, 

3, 4, 5. Then trajectory discrepancies, 

D(K) = XT(K) - XE(K) , (V-3) 

are defined for each trajectory as the difference between the final 
x-coordinate XT(K) listed as data and the final x-coordinate XE(K) 
found for the estimated trajectories in the estimated potential. The 
sum of the absolute discrepancy values 

5 

DD 1d(K)| (V-4) 

is calculated in order to obtain a quantitative measure of the error 
in the approximation of \p(x,y) by 4>(x,y). The object of each iteration 
is the reduction of DD by the addition of an ideal line quadrupole and 
its image with a magnitude and spatial location that is optimal in some 
regard. In particular, let the quadrupole be located a fixed distance 
from the region of interest (e.g., y^ = 1.4DY) and let the trajectory 
having the largest absolute discrepancy, [d(K^) 1, define the x- 
coordinate location of the quadrupole (e.g., x^ = XM(K^)). Then a quad- 
rupole magnitude A can be chosen which reduces the largest absolute 

q 

trajectory discrepancy to approximately zero. Toward this end, con- 
sider the equations of motion for electrons following the estimated 
trajectory 1C, through the estimated potential (j)(x,y). 
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8^x _ £ , 9^y _ ^ 

3t^ ■ ’ 3^2 m 3y 


(V-5) 


The integration of these equations over the time T required for elec- 
trons traveling along trajectory yields 


XE(K^) = XO(K^) + T'VXO(K^) + | | ^ 


and 


YE(Kj^) = YO(K^) + T*VY0(K^) + ^ 


3^ 

3y 


dt dt’ 


(V-6) 


(V-7) 


Recall that the trajectory consists of a number of points whose 
coordinates were determined by a stepwise solution of the equations of 
motion. Furthermore, the electric fields -|^ and were determined at 
these points by interpolation from a grid in order to proceed with the 
stepwise development of the trajectory. Since these electric fields 
and corresponding coordinate values are available, the integrals in 
equations (V-6) and (V-7) may be replaced by a pointwise summation over 
the trajectory K^. Using the trapezoidal rule for numerical integration, 
we obtain 


VXI(Kj^) = 


T 

0 


t' 


3x 


dt dt' 


N 
= Z 
n=l 


^ (S + S JAt 
2 n n-1 n 


(V-8) 


where 


S 

n 


n 
= Z 
i=l 


1 [(M) + (ii) 

2 ^3x. T 

i 1-1 




(V-9) 


Similarly, we obtain 



53 


VYI 


(K^) = 


rt' .. N 

dt dt' = Z 4 (V 


3y 


n=l ^ 


+V ^)At 
n n-1 n 


where 


(V-10) 


V 

n 


= 2 i r 

2 ^3y^ 


lAt, 


i-1 


(V-11) 


and N represents the number of time steps At^ used in the development 
of trajectory K^. Substituting these expressions for the integrals in 
equation (V-6) and (V-7) , we obtain estimates of the time T required for 
electrons following trajectory and the electron exit coordinate 
XE(K^). With YE(K^) = YO(K^) = 0, we have 

1 q 

■ 2 « VYO(Kj^) (V-12) 


and 


XE(K^) = XO(K^) + T'VXO(K^) 


+ iavxi(K^). 


(V-13) 


Notice that T > 0 for negatively charged electrons. 

We wish to add a quadrupole field <p(x,y) altering the trajectory 
only slightly and shifting the exit coordinate from its estimated 
value XE(Kj^) to the data value XT(K^). The equations of motion for 
electrons in this new potential field are 


Iji = ^ /Ji + 
3t2 m 9x 


34>, 


3x 


and 


3t2 -Sy 



(V-14) 


The integration of the first of these equations yields 


XT(K*) + X0(K*) + T'VXO(kJ) + | ^ 
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When trajectory is not too different from the trajectory, the 
* 

time T is not much different from the time T and the integral ex- 
pression in equation (V-15) can be evaluated by summation of electric 
fields over the same N points used for the evaluation of the integral 
in equation (V-6) . Therefore, 

XT(K*)=: XO(K^) + T*VX0(Kj^) + J J [VXI(K^) + Vl(K^)] (V-16) 

where 

VI(K^) ■ E i (V + V )At^ (V-17) 

n=l ^n ^-1 


and 



3(f) 3(j) 


(V-18) 


Subtracting equation (V-13) from equation (V-16) , and solving for A , 

*1 1 

we find 


^ XT(K^) - XE(K^) 

\ q ^ VI(K^) ^ 


(V-19) 



ideal line quadrupole which, when located above the trajectory ex- 
hibiting the largest discrepancy, reduces that discrepancy approximately 
to zero. 

After a quadrupole field with the appropriate magnitude is added 
to the estimated potential, a scaling operation is executed to complete 
one iteration of the scheme. The estimated potential is scaled by 
an amount sufficient to reduce the discrepancy of the scaling trajectory 
approximately to zero. To derive an algorithm for this purpose, we 



assume that electrons following the scaling trajectory experience no 
net force in the x-direction. More specifically, this implies that 
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with 


XE(6) = X0(6) + T-VX0(6) 

^ _ l£ V YI(6) 

^ 2 m VY0(6)' 


(V-20) 


(V-21) 


Then, scaling the potential by a factor SF, we shift the exit coordinate 
to a value 

XT(6) = X0(6) +T*-VX0(6) (V-22) 


where 


* 

T 


_ _ 1 i SF 

2m VY0(6) 


SF-T 


(V-23) 


Therefore, using equations (V-20) to (V-22) to solve equation (V-23) 
for the scale factor, we obtain 


_ T*_ XT(6) - X0(6) 

T ~ XE(6) - X0(6)' 


(V-24) 


Recall that this same .algorithm was used for scaling the estimated 
potential when solving the model problem by dipole synthesis. 

The main routine QP0LE.F4 implements the iterative scheme des- 
cribed above for the quadrupole synthesis of solutions to the model 
problem. QP0LE.F4 is similar to the dipole routines PTENT.F4 and 
RTENT.F4, particularly in regard to the iterative procedures used for 
scaling the estimated potential and for determining an optimal pole 
magnitude after a pole location has been specified. The major 
difference between the dipole and quadrupole routines is the criteria 
used for selecting optimal pole locations. Table 9 lists error 
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Table 9 

Error Estimates for Quadrupole Synthesis with Routine QP0LE.F4 

Average Maximum 


Iteration 

Discrepancy 

Average 

Error 

Absolute 

Error 

Absolute 

Error 


(DD) 

(%) 

(%) 

(%) 

0 

1.51 

-4.6 

5.8 

20.0 

1 

0.99 

-3.4 

4.5 

19.0 

2 

1.00 

-1.4 

4.1 

17.5 

3 

0.97 

-1.9 

4.1 

17.8 

4 

0.60 

-1.5 

2.7 

16.2 

5 

0.56 

-1.7 

2.7 

16.5 

6 

0.55 

-2.0 

2.7 

16.7 

7 

0.48 

-2.5 

3.1 

17.6 

8 

0.44 

-2.9 

3.3 

18.1 

9 

0.51 

-0.9 

3.0 

16.3 

10 

0.49 

-1.4 

3.0 

16.8 

11 

0.47 

-2.0 

3.0 

17.2 

12 

0.45 

-2.5 

3.1 

17.7 

13 

0.43 

-3.1 

3.4 

18.1 

14 

0.22 

-1.5 

1.9 

13.5 

15 

0.20 

-1.7 

2.0 

13.6 
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estimates obtained after 15 iterations of the scheme for qxiadrupole 
synthesis which can be compared to the equivalent error estimates for 
dipole synthesis appearing in Tables 4 and 5. Note that in 15 iter- 
ations of the quadrupole routine, the discrepancy is reduced to a 
value considerabley lower than the discrepancy values achieved with 15 
iterations of the dipole schemes. This discrepancy reduction for the 
quadrupole scheme is accomplished in 30’ 32" of CPU time representing 
only a small improvement in the amount of time expended per iteration. 
Although convergence limits are not reached in 15 iterations of the 
QP0LE.F4 routine, the indication is that the method would ultimately 
converge to a better approximation of the biased fin potential than 
solutions obtained with the dipole methods. Tables 10, 11 and 12 show 
the error percentages, after 5, 10 and 15 iterations respectively, 
associated with the estimated potential at 110 points chosen arbit- 
rarily from a much finer grid (i.e., the interpolation of electric fields 
with subroutine NEWTON requires a finer grid spacing) in the reference 
frame of the biased fins. Notice that the error percentages for points 
near spatial regions probed by the electron trajectories (e.g., 30<I<54 
and 5<J<17) decrease with successive iterations of the quadrupole 
scheme. The improvement of the estimated potential in these regions 
is accomplished while a degradation in the accuracy of the estimated 
potential results for points further from the space probed by the 
electron trajectories. This result for quadrupole synthesis of solu- 
tions to the model problem can be compared to similar results for dipole 
synthesis which have been discussed previously with references to Tables 
4, 5 and 6. Figure 13 shows equipotential lines for the QP0LE.F4 rou- 
tine which can be compared to Figures 8, 10 and 11 depicting the test 
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Table 10 

Estimated Potential Error Percentages for Quadrupole 
Synthesis (5 Iterations of Routine QP0LE.F4) 


Y 

21 

-16.5 

- 6.8 

0.9 

5.0 

5.5 

2.9 

0.4 

1.1 

-0.4 

- 1.6 

- 5.7 

C 

19 
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- 0.6 
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0 

0 
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0.5 

-0.2 
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R 

15 
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0.8 
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D 
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N 

D 
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Table 11 

Estimated Potential Error Percentages for Quadrupole 
Synthesis (10 Iterations of Routine QP0LE.F4) 


Y 
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-16.8 

- 7.1 

0.1 3.5 
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0.7 
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2.2 
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0.9 

0.5 
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D 

I 

13 
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0.9 
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N 

11 
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1.4 
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1.0 

1.8 
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T 
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0.9 
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Table 12 

Estimated Potential Error Percentages for Quadrupole 
Synthesis (15 Iterations of Routine QP0LE.F4) 


Y 

21 

-12.4 

- 1.6 

2.7 

-0.6 

-1.4 

-0.1 

0.3 
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0.4 
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C 
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0.5 
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0.1 
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potential and the dipole synthesized solutions. 

B. Potential Synthesis with Quadrupole Arrays 

In order to speed the process of discrepancy reduction,, consider 
the placement of five ideal line quadrupoles upon each iteration of a 
new scheme for quadrupole synthesis of solutions to the model problem. 
The procedure for calculating an initial estimated potential, 4>^(x,y), 
based upon a uniform field approximation and the data for the scaling 
trajectory (i.e., K=6) is unchanged and described by equations IV-1) 
through (IV-4) . The estimated potential can then be modified by an 
iterative scheme that locates five quadrupoles simultaneously in order 
to achieve discrepancy reduction for all five of the electron tra- 
jectories. Equation (V-25) depicts this situation 

N 5 

<p(x,y) = <{) (x,y) + E E (j)^(x,y) (V-25) 

° 1=1 K =1 ^ 

where <f) (x,y) is given by equation (V-1) and the i index runs up to 
K 

the total number of iterations that are needed for adequate discrepancy 
reduction. In developing this scheme, we expect that each Iteration 
will require more CPU time since five quadrupoles are located as com- 
pared to only one with the old scheme. However, we hope that fewer 
iterations will be required to achieve a given reduction in discrepancy 
and that the total CPU time expended will decrease. To define the new 
scheme, let all five quadrupoles be located with a y-coordinate value 
of 1.4DY. Let the x-coordinate values of the quadrupole locations be 
equal to the x-coordinates corresponding to the highest points on the 
five estimated trajectories that have been developed for the estimated 


potential. 
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= XM(K) ; 1,2, 3. 4, 5 


(V-26) 


Once the five quadrupole locations have been chosen, we seek five 
quadrupole Biagnitudes (i.e. , A^;K = 1, 2, 3, 4, 5) which reduce the 
discrepancies associated with all five trajectories to approximately 
zero. Toward this end, consider the equations of motion for electrons 
moving in the estimated potential <j)(x,y). Then, in analogy with equa- 
tions (V-5) through (V-13) , an equation for each estimated trajectory 
can be determined. 


XE(K) = X0(K) + T.VXO(K) + ^ VXI(K) ; K = 1,2, 3, 4, 5 (V-27) 

2 m 

Now, we wish to add five quadrupole fields to the estimated potential 
so as to alter the five estimated trajectories only slightly and shift- 
ing the estimated exit coordinates (XE(K) ; K=l, 2, 3, 4, 5) over to 
the exit coordinates (XT(K) ; K= 1, 2, 3, 4, 5) measured as data. 

The equations of motion in this new estimated potential are 


and 



3i 


5 

+ L 
K=1 


3(J) 


K. 


3x 


(V-28) 


3t2 m ^9y 


5 

+ Z 
K=1 




K 


3x 


). 


(V-29) 


An integration of the equations of motion subject to the initial and 
final conditions available as trajectory data results in five equations 
describing the electron trajectories in a new estimated potential. 

Thus, 

, 5 

XT(K) = X0(K) + T-VXO(K) +4-^ [VXI(K) + Z 

2 m , 


A^VI(K)], (V-30) 


with VI (K) defined by equations (V-17) and (V-18) and where 



K = 1, 2, 3, 4, 5. The subtraction of equation (V-27) from equation 
(V-30) yields. 
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XT(K) - XE(K) =1^ A^VI(K) ; K = 1,2,3, 4, 5. (V-31) 

Equation (V-31) represents a linear system of five equations in 
five unknowns (i.e. , = 1, 2, 3, 4, 5) which can be solved by 

gaussian elimination procedures to determine the magnitude of the five 
ideal line quadrupoles. When the potential fields associated with 
these quadrupoles are added to the estimated potential, a new estimated 
potential is synthesized ^ Then the discrepancies associated with the 
new estimated trajectories, developed in the new estimated potential, 
should individually be approximately equal to zero. Thus, complete 
discrepancy reduction could be accomplished with a single application 
of the new scheme. However, the scheme involves several approximations 
which produce a less than total reduction of discrepancy. For example, 
the time T in equations (V-27) and (V-30), as well as the VXI(K) 
terms, depend upon the old and new estimated potentials respectively. 
Therefore, the terms involving these expressions will not cancel ex- 
actly as indicated in the derivation of equation (V-31) . The approx- 
imation as implied by equations (V-27) and (V-30) is that the five 
electron trajectories are not altered greatly by the addition of five 
quadrupoles for the reduction of discrepancy. Therefore, the errors 
introduced by these approximations will result in a less than total 
reduction of discrepancy for any one iteration of the scheme. However, 
on successive iterations, the approximations will produce lesser errors 
and for a suitable initial potential estimate (J)^, the method of 
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successive approximation will result in a convergent iterative scheme. 

The main routine QP0LES.F4 (Appendix E) implements the method 
described above for the synthesis of solutions to the model problem. 

Table 13 shows the error estimates generated at 110 points on a grid 
located between the biased fins and above the zero potential reference 
plane. These results can be compared to Tables 10, 11 and 12 which 
list the results of routine QP0LE.F4 also used for quadrupole synthesis 
of solutions to the model problem. Note that the single pole scheme 
results in solutions that are scaled better than the solutions syn- 
thesized by the addition of five poles simultaneously. However, 
discrepancy reduction is much more rapid for the latter method. A 
reduction in discrepancy from 1.5 to 0.4 can be accomplished with a 
single iteration of the QP0LES.F4 scheme requiring 3 '24" of computer 
processing time (CPU time) as compared to 14 iterations of the QP0LE.F4 
scheme requiring 28*29". Although 4 iterations of the method have been 
listed, the reduction of discrepancy and of approximation errors reaches 
a limit after only 2 iterations. This is in contrast to the single pole 
scheme where convergence limits had not been reached after 15 iterations. 
Figure 14 shows equipotential lines for the QP0LES.F4 routine. 

Although convergence is rapid for the main routine QP0LES.F4, the 
synthesized solution to the model problem underestimates the biased fin 
potential at 109 of the 110 points considered. The scaling procedure 
itself is unchanged from the successful procedure used previously for 
the single quadrupole and dipole routines. The difficulty seems to be 
manifested in balancing the strategies of field scaling and field shap- 
ing depending on whether the addition of quadrupoles, reducing discrep- 
ancies for the first five trajectories, has a sufficiently detrimental 
effect on the discrepancy of the scaling trajectory. To possibly save 
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Table 13 


Error Estimates for Routine QP0LES.F4 with Step Size d 
(0.013 <d< 0.025) and Scale Tolerance of 1.0% 


Iterations 

Discrepancy 

Scale 

Error 

Average 

Error 

Average 

Absolute 

Error 

Maximum 

Absolute 

Error 


(DD) 

(%) 

(%) 

(%) 

(%) 

0 

1.51 

0.05 

-4.6 

5.8 

20.0 

1 

0.38 

0.02 

-3.2 

3.3 

14.0 

2 

0.08 

0.16 

-3.3 

3.3 

13.7 

3 

0.05 

0.02 

-3.2 

3.3 

14.1 

4 

0.06 

0.22 

-3.3 

3.3 

14.5 
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17 - 5.0 - 3.3 - 1.5 -1.4 -1.5 -1.6 -1.6 -1.9 - 2.5 - 3.3 - 3.0 

15 - 3.8 - 3.2 - 2.0 -1.7 -1.7 -1.7 -1.8 -2.0 - 2.5 - 2.8 - 2.0 

13 - 3.1 - 3.0 - 2.3 -1.9 -1.8 -1.9 -2.1 -2.4 - 2.4 - 2.4 - 1.4 

11 - 2.8 - 3.0 - 2.5 -2.1 -1.9 -1.9 -2.0 -2.2 - 2.4 - 2.3 - 1.3 

9 - 2.8 - 3.1 - 2.8 -2.4 -2.1 -2.0 -2.2 -2.4 - 2.6 - 2.4 - 1.4 

7 - 3.4 - 3.6 - 3.3 -2.8 -2.3 -2.2 -2.4 -2.8 - 3.1 - 2.9 - 2.0 

5 - 5.1 - 5.4 - 4.7 -3.6 -2.7 -2.4 -2.7 -3.5 - 4.3 - 4.6 - 3.9 

3 -14.2 -13.7 -11.0 -7.2 -3.9 -2.6 -3.9 -7.2 -10.7 -13.0 -13.0 
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time, a decision to scale or not is made by comparing the scale error 
(i.e., discrepancy) of this trajectory with an arbitrary scale tolerance. 
Tables 14, 13 and 15 show error estimates for the synthesis of 
solutions to the model problem for scale tolerances of 5.0%, 1.0% and 
0.1% respectively. Note that a scale tolerance of 0.1% is preferable. 
Decreasing the scale tolerance further will not Improve the accuracy 
of the solution, since the routine already scales at every opportunity 
with a scale tolerance of 0.1%. Thus, invoking the scaling operation 
after each shaping operation facilitates the reduction of discrepancy 
and produces a better approximation of the biased fin potential. How- 
ever, the synthesized solution to the model problem still underestimates 
the biased fin potential at 109 of the 110 points considered. Thus, 
field shaping by the addition of five quadrupoles simultaneously re- 
sults in a rapid convergence at the expense of accurate scaling. 

Recall that the terms VI (K) and VXI(K) in equations (V-17) and 
(V-8) are calculated by summing terms generated by the trapezoidal 
rule for integration over the estimated electron trajectories. The 
accuracy of this integration depends upon the number of trapezoidal 
terms as determined by the size of the time step used for tracking 
electrons through a potential field. Subroutine DVOGEL monitors the 
trajectory arc length corresponding to a particular time step to 
determine the acceptance or rejection of that time step. Whenever the 
arc length surpasses some upper limit, the time step is reduced and 
conversely whenever the arc length falls below some lower limit, the 
time step is increased. Table 16 shows x and y coordinate values at 
points along a particular electron trajectory as well as the electric 
field at those points for an ideal line quadrupole that will be added 
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Table 14 

Error Estimates for Routine QP0LES.F4 with Scale Tolerance of 5.0% 


Iterations 

Discrepancy 

Scale 

Error 

' 

(DD) 

(%) 

0 

1.51 

0.05 

1 

0.62 

2.78 

2 

0.26 

2.64 

3 

0.05 

2.51 

4 

0.03 

2.51 


Average Maximum 

Average Absolute Absolute 

Error Error Error 

(%) (%) (%) 

-4.6 5.8 20.0 

-7.0 7.0 17.4 

-6.9 6.9 17.3 

-6.8 6.8 17.3 

-6.7 6.7 17.6 
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T 
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D 

E 

X 
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21 

17 

13 
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21 -16.2 - 6.5 

19 -11.3 - 7.4 

17 - 9.0 - 7.3 

15 - 7,8 - 7.1 

13 - 7.1 - 6.8 

11 - 6.7 - 6.7 

9 - 6.7 - 6.9 
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Table 15 


Error Estimates for Routine QP0LES.F4 with Scale Tolerance of 0.1% 


Iterations 

Discrepancy 

Scale 

Error 

Average 

Error 

Average 

Absolute 

Error 

Maximum 

Absolute 

Error 


(DD) 

(%) 

(%) 

(%) 

(%) 
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1.51 

0.05 

-4.6 

5.8 
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0.38 
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STEP 


Table 16 

Quadrupole Electric Fields at Points Defining Estimated 
Electron Trajectories in the Uniform Field Potential 


STEP SIZE 


.008 

<d<.025 

.013 

<d< .025 

.013 

<d < .050 

Y 

E 

Y 

E 

Y 

E 


X 


X 


X 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.03 

-0.05 

0.05 

-0.08 

0.05 

-0.08 

0.07 

-0.11 

0.10 

-0.16 

0.10 

-0.16 

0.10 

-0.16 

0.15 

-0.24 

0.15 

-0.24 

0.13 

-0.21 

0.20 
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0.20 

-0.33 

0.16 
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0.24 

-0.42 

0.24 

-0.42 

0.20 

-0.33 

0.29 
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0.29 
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0.23 

-0.40 

0.34 

-0.66 

0.34 

-0.66 

0.26 

-0.46 

0.39 

-0.80 

0.39 

-0.80 

0.30 

-0.54 

0.44 

-0.97 

0.44 

-0.97 

0.33 

-0.62 

0.49 

-1.17 

0.49 

-1.17 

0.36 

-0.71 

0.54 

-1.41 

0.54 

-1.41 

0.39 

-0.81 

0.59 

-1.73 

0.59 

-1.73 

0.42 

-0.92 

0.64 

-2.13 

0.64 

-2.13 

0.46 

-1.05 

0.69 

-2.65 

0.69 

-2.65 

0.49 

-1.19 

0.74 

-3.33 

0.74 

-3.33 

0.53 

-1.35 

0.78 

-4.18 

0.78 

-4.18 

0.56 

-1.54 

0.83 

-5.43 

0.83 

-5.43 

0.59 

-1.76 

0.88 

-7.13 

0.88 

-7.13 

0.62 

-2.01 

0.82 

-5.03 

0.82 

-5.03 

0.65 

-2.31 

0.73 

-3.30 

0.65 

-2.24 

0.69 

-2.68 

0.64 

-2.18 

0.48 

-1.13 

0.72 

-3.08 

0.55 

-1.50 

0.30 

-0.56 

0.75 

-3.57 

0.46 

-1.06 

0.12 

-0.19 

0.78 

-4.17 

0.38 

-0.75 



0.82 

-4.94 

0.29 

-0.53 



0.85 

-5.89 

0.20 

-0.33 



0.88 

-7.04 

0.11 

-0.17 



0.84 

-5.73 

0.02 

-0.04 



0.75 

-3.57 





0.66 

-2.38 





0.57 

-1.65 





0.48 

-1.14 





0.40 

-0.82 





0.30 

-0.56 





0.21 

-0.36 





0.12 

-0.20 





0.04 

-0.06 
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for field shaping. This trajectory shows no variation in the x coor- 
dinate since the estimated potential field used to generate the tra- 
jectory from the initial electron data is the uniform field approximation 
developed as the first operation in all of the iterative routines that 
have been discussed. Note that decreasing the lower limit on the step 
size results in more steps as the electron moves against the potential 
gradient (increasing y coordinate) . Note also that increasing the upper 
limit on the step size results in fewer steps as the electron moves with 
the potential gradient (decreasing y coordinate) . Tables 17 and 18 when 
compared to Table 15 show that the three different bounds on the step 
size result in synthesized solutions to the model problem which are 
only slightly different from each other. 

It is interesting to consider the effects of erroneous measurements 
of initial electron energy or initial electron direction. Tables 19 and 
20 show error estimates for the approximations of the biased fin poten- 
tial which result from erroneous measurements of electron initial 
energies of 3% and 6% respectively. Note that this type of error affects 
the field scaling operation but has little effect on the field shaping 
operations. Tables 21 and 22 show error estimates for approximations 
of the biased fin potential resulting from erroneous measurements of 
electron initial directions of 1® and 3“ respectively. Note that this 
type of error severely affects field shaping operations as well as the 


scaling operations. 
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Table 17 




Error Estimates 

for Routine 

QP0LES.F4 with Step Size 

d 






(0.008 <d< 

0.025) 












Average 

Maximum 







Scale 

Average 

Absolute 

Absolute 

Iterations 

Discrepancy 

Error 


Error 

Error 

Error 





(DD) 


(%) 


(%) 

(%) 

(%) 


0 



1.51 


0.05 


-4.6 

5.8 

20.0 


1 



0.38 


0.02 


-3.2 

3.3 

14.1 


2 



0.09 


0.22 


-3.4 

3.5 

13.9 


3 



0.04 


0.18 


-3.3 

3.4 

14.1 


4 



0.02 


0.21 


-3.3 

3.4 

14.3 

Y 


21 

-12.7 

- 2.1 

1.5 

-1.1 - 

•1.4 

-1.3 -1.3 

-1.1 - 2.6 

-3.8-8 

C 


19 

- 7.6 

- 3.3 

- 0.5 

-1.0 - 

•1.4 

-1.4 -1.4 

-1.6 - 2.6 

-3.8-4 

0 

0 


17 

- 5.3 

- 3.4 

- 1.5 

-1.3 - 

•1.4 

-1.5 -1.6 

-1.9 - 2.6 

-3.3-3 

R 


15 

- 4.0 

- 3.2 

- 2.0 

-1.6 - 

•1.6 

-1.6 -1.7 

-2.0 - 2.5 

-2.8-2 

D 

I 


13 

- 3.3 

- 3.1 

- 2.3 

-1.8 - 

•1.7 

-1.7 -1.9 

-2.1 - 2.5 

-2.5-1 

N 


11 

- 3.0 

- 3.1 

- 2.5 

-2.1 - 

•1.9 

-1.9 -2.0 

-2.2 - 2.5 

-2.3-1 

A 

T 


9 

- 3.0 

- 3.2 

- 2.8 

-2.4 - 

•2.1 

-2.0 -2.1 

-2.4 - 2.6 

-2.4-1 

E 


7 

- 3.6 

- 3.8 

- 3.4 

-2.8 - 

•2.3 

-2.2 -2.4 

-2.8 - 3.1 

-3.0-2 

I 


5 

- 5.3 

- 5.5 

- 4.7 

-3.6 - 

•2.7 

-2.3 -2.7 

-3.5 - 4.4 

-4.7-3 

N 

D 


3 

-14.4 

-13.8 

-11.0 

-7.2 - 

•3.9 

-2.5 -3.9 

-7.2 -10.7 

-13.0 -13 


E 

X 

(J) 


22 26 30 34 38 42 46 

X COORDINATE INDEX (I) 


50 


54 


58 


62 
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Table 18 


Error Estimates for Routine QP0LES.F4 with Step Size d 

(0.013 <d< 0.050) 


Iterations 

Discrepancy 

Scale 

Error 

Average 

Error 

Average 

Absolute 

Error 

Maximum 

Absolute 

Error 


(DD) 

(%) 

(%) 

(%) 

(%) 

0 

1.51 

0.43 

-6.9 

6.9 

17.4 

1 

1.01 

0.01 

-3.2 

3.3 

14.2 

2 

0.85 

0.02 

-3.9 

4.3 

17.3 

3 

0.28 

0.03 

-3.1 

3.3 

14.8 

4 

0.08 

0.37 

-3.3 

3.4 

15.1 


Y 

21 

-13.7 

- 3.3 

0.2 

-2.0 

-1.9 

-0.6 

-0.5 

1.1 

- 1.9 

- 4.2 

- 7.8 

C 

19 

- 8.6 

- 4.4 

- 1.5 

-1.8 

-1.8 

-0.9 

-0.6 

-0.2 

- 1.9 

- 3.6 

- 3.9 

0 

0 

17 

- 6.3 

- 4.4 

- 2.4 

-2.0 

-1.8 

-1.2 

-0.9 

-0.9 

- 1.8 

- 2.8 

- 2.2 

R 

15 

- 5.0 

- 4.2 

- 2.9 

-2.2 

-1.9 

-1.4 

-1.2 

-1.2 

- 1.8 

- 2.2 

- 1.3 

D 

I 

13 

- 4.2 

- 4.0 

- 3.1 

-2.4 

-2.0 

-1.6 

-1.4 

-1.4 

- 1.7 

- 1.8 

- 0.8 

N 

11 

- 3.9 

- 3.9 

- 3.3 

-2.6 

-2.1 

-1.8 

-1.6 

-1.6 

- 1.8 

- 1.6 

- 0.6 

A 

T 

9 

- 3.9 

- 4.1 

- 3.6 

-2.9 

-2.3 

-1.9 

-1.8 

-1.8 

- 2.0 

- 1.7 

- 0.7 

E 

7 

- 4.4 

- 4.6 

- 4.1 

-3.3 

-2.5 

-2.1 

-2.0 

CM 

• 

CM 

1 

- 2.4 

- 2.3 

- 1.8 

I 

5 

- 6.2 

- 6.3 

- 5.4 

-4.1 

-2.9 

-2.3 

-2.3 

-3.0 

- 3.7 

- 3.9 

- 3.2 

N 

D 

T? 

3 

-15.1 

-14.5 

-11.7 

-7.7 

-4.1 

-2.5 

-3.6 

-6.6 

-10.1 

-12.4 

-12.4 

X 


22 

26 

30 

34 

38 

42 

46 

50 

54 

58 

62 


X COORDINATE INDEX (I) 
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Table 19 


Error Estimates for Routine QP0LES.F4 with 3% Error in Beam Energy 


Iterations 

Discrepancy 

Scale 

Error 

Average 

Error 

Average 

Absolute 

Error 

Maximum 

Absolute 

Error 


(DD) 

(%) 

(%) 

(%) 

(%) 

0 

1.51 

0.05 

-1.7 

5.0 

17.6 

1 

0.36 

0.02 

-0.1 

1.8 

10.9 

2 

0.09 

0.18 

-0.2 

1.6 

10.7 

3 

0.03 

0.17 

-0.2 

1.6 

11.0 

4 

0.01 

0.18 

-0.2 

1.6 

11.1 


Y 

21 

- 8.9 

2.7 

5.7 

1.4 

0.5 

0.7 

0.5 

1.0 

0.2 

0.0 

- 4.7 

C 

19 

- 3.6 

1.0 

3.4 

1.8 

0.9 

0.7 

0.6 

0.7 

0.4 

- 0.2 

- 0.7 

0 

0 

17 

- 1.3 

0.6 

2.2 

1.7 

1.0 

0.7 

0.7 

0.6 

0.4 

0.3 

0.9 

R 

15 

- 0.1 

0.6 

1.6 

1.5 

1.0 

0.8 

0.7 

0.6 

0.6 

0.7 

1.8 

D 

I 

13 

0.5 

0.7 

1.2 

1.2 

1.0 

0.8 

0.7 

0.6 

0.6 

1.0 

2.2 

N 

11 

0.8 

0.6 

0.9 

1.0 

0.9 

0.8 

0.6 

0.6 

0.6 

1.1 

2.4 

A 

T 

9 

0.7 

0.4 

0.6 

0.7 

0.8 

0.7 

0.5 

0.4 

0.5 

1.0 

2.2 

E 

7 

0.1 

- 0.3 

-2.7 

0.3 

0.6 

0.6 

0.4 

0.1 

0.0 

0.4 

1.6 

I 

5 

- 1.8 

- 2.0 

-1.4 

-0.5 

0.2 

0.4 

0.1 

• 

0 

1 

-1.3 

- 1.3 

- 0.3 

N 

D 

T? 

3 

-11.1 

-10.6 

-8.0 

-4.3 

-1.0 

0.2 

-1.2 

-4.4 

-7.8 

-10.0 

-10.0 

X 


22 

26 

30 

34 

38 

42 

46 

50 

52 

58 

62 


X COORDINATE INDEX (I) 



30 


54 
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Table 20 


Error Estimates for Routine QPOLES.F4 with 6% Error in Beam Energy 


Iterations 

Discrepancy 

Scale 

Error 

Average 

Error 

Average 

Absolute 

Error 

Maximum 

Absolute 

Error 


(DD) 

(%) 

(%) 

(%) 

(%) 

0 

1.51 

0.05 

1.2 

5.2 

15.1 

1 

0.38 

0.02 

2.7 

3.8 

8.8 

2 

0.07 

0.16 

2.6 

3.6 

8.5 

3 

0.04 

0.14 

2.7 

3.7 

8.8 

4 

0.02 

0.15 

2.7 

3.7 

9.0 


Y 

21 

-7.2 

4.1 

7.8 

4.9 

4.6 

4.7 

4.9 

5.1 

3.6 

2.3 

-3.2 

C 

19 

-1.7 

2.8 

5.7 

5.0 

4.7 

4.7 

4.7 

4.5 

3.5 

2.2 

1.2 

0 

0 

17 

0.7 

2.6 

4.6 

4.8 

4.6 

4.6 

4.5 

4.2 

3.5 

2.7 

3.0 

R 

15 

2.0 

2.8 

4.0 

4.5 

4.5 

4.5 

4.4 

4.1 

3.6 

3.3 

4.0 

D 

I 

13 

2.7 

3.0 

3.8 

4.2 

4.4 

4.4 

4.2 

4.0 

3.6 

3.6 

4.6 

N 

11 

3.1 

3.0 

3.5 

4.0 

4.2 

4.2 

4.1 

3.8 

3.6 

3.8 

4.8 

A 

T 

9 

3.0 

2.8 

3.2 

3.7 

4.0 

4.1 

3.9 

3.7 

3.4 

3.6 

4.7 

E 

7 

2.5 

2.2 

2.6 

3.3 

3.7 

3.9 

3.7 

3.3 

3.0 

3.1 

4.0 

I 

5 

0.6 

0.4 

1.2 

2.4 

3.1 

3.7 

3.1 

2.5 

1.6 

1.3 

2.9 

N 

D 

E 

3 

-9.0 

-8.9 

-5.5 

-1.5 

2.1 

3.5 

2.1 

-1.4 

-5.1 

-7.6 

-7.7 

X 

;j) 


22 

26 

30 34 38 42 

X COORDINATE INDEX 

46 

(I) 

50 

54 

58 

62 



30 


42 


54 
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Table 21 


Error Estimates for Routine QP0LES.F4 



with 1% Error 

in Initial 

Electron 

Direction 


Iterations 

Discrepancy 

Scale 

Error 

Average 

Error 

Average 

Absolute 

Error 

Maximum 

Absolute 

Error 


(DD) 

(%) 

(%) 

(%) 

(%) 

0 

1.57 

0.05 

-4.9 

6.0 

20.3 

1 

0.43 

0.01 

-3.7 

3.8 

-15.0 

2 

0.11 

0.37 

-4.0 

4.0 

-15.7 

3 

0.06 

0.38 

-3.9 

3.9 

-16.2 

4 

0.03 

0.44 

-4.0 

4.0 

-16.5 


4.7 0.1 -1.5 -1.9 -1.7 -1.9 -1.4 - 2.9 - 3.4 

5.8 - 2.0 -1.8 -2.0 -1.9 -1.9 -1.9 - 2.6 - 3.1 

5.8 - 3.1 -2.3 -2.1 -2.0 -2.0 -2.0 - 2.4 - 2.4 

5.6 - 3.7 -2.7 -2.3 -2.1 -2.1 -2.1 - 2.2 - 1.9 

5.4 - 4.0 -3.0 -2.5 -2.3 -2.2 -2.1 - 2.1 - 1.6 

5.3 - 4.3 -3.3 -2.7 -2.4 -2.3 -2.2 - 2.1 - 1.5 

5.4 - 4.6 -3.6 -3.0 -2.6 -2.4 -2.4 - 2.2 - 1.6 

5.9 - 5.1 -4.1 -3.2 -2.7 -2.6 -2.7 - 2.7 - 2.2 

7.5 - 6.4 -4.9 -3.6 -2.9 -2.9 -3.5 - 3.9 - 3.9 

15.6 -12.6 -8.5 -4.8 -3.1 -4.1 -7.1 -10.3 -12.3 

26 30 34 38 42 46 50 54 58 

X COORDINATE INDEX (I) 





30 


42 


54 
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Table 22 

Error Estimates for Routine QP0LES.F4 



with 3% Error 

in Initial 

Electron 

Direction 


Iterations 

Discrepancy 

Scale 

Error 

Average 

Error 

Average 

Absolute 

Error 

Maximum 

Absolute 

Error 


(DD) 

(%) 

(%) 

(%) 

(%) 

0 

1.83 

0.05 

-6.0 

6.7 

21.2 

1 

0.58 

0.04 

-4.9 

5.0 

20.1 

2 

0.21 

0.80 

-5.6 

5.6 

21.9 

3 

0.13 

0.05 

-4.3 

4.7 

21.6 

4 

0.07 

0.31 

-4.5 

5.1 

-22.3 


21 -22.3 - 9.6 - 2.3 - 1.6 -2.4 -2.6 -2.7 -2.6 -2.9 -2.1 -0.1 

19 -16.2 -10.3 - 4.6 - 2.8 -2.7 -2.7 -2.6 -2.4 -2.1 -0.8 2.7 

17 -13.1 -10.0 - 5.8 - 3.7 -3.1 -2.8 -2.5 -2.1 -1.4 0.2 3.5 

15 -11.3 - 9.4 - 6.5 - 4.4 -3.4 -2.9 -2.5 -1.9 -1.0 0.8 3.8 

13 -10.1 - 9.0 - 6.8 - 4.9 -3.7 -3.0 -2.4 -1.7 -0.6 1.2 3.9 

11 - 9.5 - 8.7 - 7.0 - 5.3 -4.0 -3.1 -2.5 -1.7 -0.5 1.3 3.8 

9 - 9.2 - 8.8 - 7.3 - 5.6 -4.3 -3.3 -2.5 -1.7 -0.6 1.1 3.4 

7 - 9.6 - 9.2 - 7.8 - 6.0 -4.5 -3.4 -2.7 -2.0 -1.0 0.5 2.6 

5 -11.1 -10.7 - 9.0 - 6.9 -4.9 -3.6 -3.0 -2.7 -2.3 -1.2 0.6 

3 -19.5 -18.5 -15.1 -10.4 -6.1 -3.8 -4.2 -6.3 -8.7 -9.9 -9.0 

22 26 30 34 38 42 46 50 54 58 62 

X COORDINATE INDEX (I) 


21 

17 
13 

9 
5 

30 42 54 



Y 

C 

0 

0 

R 

D 

I 

N 

A 

T 

E 

I 

N 

D 

E 

X 

(J) 
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VI. EXPERIMENTAL RESULTS 

The biased fin potential referred to in sections III, IV and V of 
this report has been a conceptual entity rather than a physical entity. 
The same can be said for the electron trajectories that are used to 
measure and probe the biased fin potential. In this section, however, 
we will describe a laboratory experiment that was conducted to test 
the computer codes that have been developed for the synthesis of 
approximations to the biased fin potential. 

The potential developed in section III was proposed to approximate 
the field existing near an infinite number of parallel-planar fins, all 
of which extended infinitely in both parallel directions. Thus, the 
conceptual system extended infinitely in three spatial directions. How- 
ever, the physical system was constructed inside a cylindrical bell jar 
and consisted of six parallel plates approximately five inches high 
and were biased with -100 volts at a distance of 1.7 inches from a zero 
potential ground plane. Due to the finite size of the system, we ex- 
pect that the physical potential field existing inside the evacuated 
bell jar will be somewhat different from the biased fin potential as 
derived in section III. In spite of this dimensional limitation, we 
hoped to use the computer codes described in section V to process 
electron trajectory data and develop a synthesized potential field for 
comparision with the biased fin potential plotted in Figure 8. 

An electron beam was generated by an electron gun consisting of a 
wire filament inside a cylindrical aluminum can having a pinhole at the 
center of one end. The electron gun was positioned below a grounded 
wire mesh (i.e., the reference plane) so that the electron beam emanating 
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from the can. would pass through the wire mesh on a traj ectory that 
ultimately returned to the wire mesh due to the deflection of the 
electron beam by the negatively charged fins. The source and sink 
points in the reference plane were recorded for six trajectories as 
listed in Table 23. Although the dimensions of the physical system 
(e.g., 3.17 inches by 1.7 inches) are different from those of the con- 
ceptual system (e.g., 1.8 units by 1 unit), it is apparent that the 
trajectories listed as physical data do not correspond to those tra- 
jectories defined conceptually in the statement of the model problem. 
Thus, it is interesting to speculate what the electron trajectories 
might have looked like in the conceptual biased fin potential given the 
initial electron parameters measured as physical data. Figure 15 shows 
projected trajectories developed for an assumption of the biased fin 
potential as expressed in equation (III-IO) and for the initial electron 
velocities listed in Table 23. Note that the sink points along the 
reference plane do not correspond to the sink points measured as phy- 
sical data. Therefore, we conclude that the physical potential field 
of the biased fins must be somewhat different from the conceptualized 
potential. With the knowledge that the unknown potential field in the 
bell jar is similar to, but not identical to, the conceptualized biased 
fin potential; it is appropriate to apply the computerized methods des- 
cribed in section V in order to synthesize approximations of the unknown 
potential. Figure 16 shows the results obtained with five iterations 
of the QP0LE.F4 routine using all the data that appears in Table 23. 

Note that this solution indicated a lower potential midway between the 
fins when compared to the potential field of Figure 8. This may well 
be the result of approximating an infinite fin height (measured 
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(I) 
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Approximation of the Experimental Potential Synthesized 
with QP0LE.F4 


J=22 



Approximation of the Experimental Potential S 3 mthesized 
with QP0LES.F4 
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perpendicular to the reference plane) with a fin only 5 inches high. 

Figure 17 shows the results obtained when only four of the tra- 
jectories (i.e. , K = 2,3,4 and 6) are used with the QP0LES.F4 routine. 
Attempts to use more of the data in Table 23 caused instabilities in 
the method preventing the synthesis of approximations to the potential 
field existing in the bell jar. Thus, by concentrating only on the 
interior trajectories, a potential field can be created which reproduces 
the sink points of those interior trajectories. However, as seen in 
Figure 17, the sink points for the K=1 and K=5 trajectories are in 
error since the data for these trajectories is not used. 

Thus, we have demonstrated that the quadrupole methods developed 
for an idealized biased fin potential and a particular set of trajectory 
data can be used to synthesize approximations of the similar potential 
field existing in the bell jar by using a different set of trajectory 
data. We see that the choice of data for a particular scheme does 
affect the convergence properties of that scheme, to the extent that no 
solution is a possibility. Therefore, we conclude that data selection 
is important for the application of these iterative routines. Con- 
versely, the definition of a particular iterative scheme may depend 
heavily on the data that can be obtained. 
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Table 23 

Experimental Electron Trajectory Data 


Traj ectory 

Source 

Points 

Source Velocities 

Sink 

Points 

K 

XO(K) 

Y0(K) 

VXO(K) 

VYO(K) 

XT(K) 

YT(K) 


(Inches) 

(inches) 

(— ) 
sec 

(— ) 
sec 

(inches) 

(inches) 

1 

0.73 

0.0 

0.0 

5.5x10^ 

2.72 

0.0 

2 

1.25 

0.0 

0.0 

5.5x10^ 

2.29 

0.0 

3 

2.27 

0.0 

0.0 

5.5x10^ 

1.52 

0.0 

4 

2.54 

0.0 

0.0 

5.5x10^ 

0.96 

0.0 

5 

2.91 

0.0 

0.0 

5.5x10^ 

0.33 

0.0 

6 

0.12 

0.0 

3.4x10^ 

3.4x10^ 

3.29 

0.0 
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VII. CONCLUSIONS 


The synthesis of electric potentials from the fields of electric 
dipoles and quadmpoles has been considered as a method for approxima- 
ting unknown potential fields in two dimensions. Although we were most 
interested in determining potential fields on an unbounded spatial do- 
main; computational methods were developed first for the classical 
boundary value problem in electrostatics, since other well-known meth- 
ods exist for its solution. 

Several iterative schemes for the dipole synthesis of solutions 
to a square boundary value problem were developed and their accuracy 
was compared to a series expansion solution. We found that as the di- 
pole singularities were moved farther away from the boundary, the syn- 
thesized solutions became more accurate at the expense of more compu- 
tation (i.e. , more iterations of a scheme were required for convergence 
to a desired limit) . Additionally, procedures for scaling, shifting 
and rotating the estimated potential were developed to assist in the 
solution of the boundary value problem. The rotational procedure was 
found to be of assistance only in determining an initial estimated 
potential that best matched the known boundary values. The scaling and 
shifting procedures were found to be helpful immediately after each 
successive addition of a dipole field. The dipole fields were chosen 
to reduce the discrepancy between a known potential value and its 
estimated value at the point on the boundary where the discrepancy was 
largest. This scheme showed good convergence properties. We found that 
convergence could be hastened by considering potential discrepancy 
reduction at two boundary points simultaneously and then allowing the 
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location of dipole field singularities at variable distances from the 
boundary. Extensions of this work on the boundary value problem could 
be interesting if quadrupole synthesis is explored and if other bound- 
ary shapes are considered. 

Both dipole and quadrupole methods have been developed for the 
synthesis of unknown potential fields in two dimensions of an unbounded 
spatial domain. Our first step was to introduce the concept of a ref- 
erence plane existing near the surface charge distribution responsible 
for the potential field existing in some spatial region of interest. 

The reference plane represents an artificial boundary condition, since 
we have required a potential of zero on it in order to allow the pres- 
ence of a measuring device (i.e. , an electron beam probe) in the vicinity 
of the charged surface. Charged particles emanate from the reference 
plane, they are deflected by the potential field of a charge distribu- 
tion and fall back to the reference plane. The energy of the charged 
particle and it's initial and final coordinate values have been used 
as constraints upon acceptable dipole and quadrupole synthesized 
approximations to the potential resulting from the charge distribution. 
Dipole and quadrupole fields have been added as image pairs to insure 
a zero potential on the reference plane. 

Basic schemes have been developed for the dipole synthesis of 
electric potentials as constrained by electron trajectory data. A 
uniform field approximation is determined so as to meet the constraining 
data for a particular trajectory called the scaling trajectory. Then 
a dipole field is added in order to meet the constraints of two other 
trajectories while keeping the dipole field singularity at a fixed 
distance from the reference plane. A potential shifting operation has 
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not been required since the dipole fields have been added in a manner 
that kept the potential at zero on the reference plane. However, 
potential scaling has been found to be beneficial after the addition 
of each dipole field in order to meet the constraints of the scaling 
trajectory. Convergence limits for this scheme have been adequate, 
although for the model problem, the synthesized solution overestimated 
the test potential. We found that an underestimation of the test 
potential could be obtained by altering the criterion used for deter- 
mining which two sets of trajectory constraints are satisfied when 
adding a dipole field. Therefore, we were able to establish bounds on 
the unknovm. potential created by the surface charge distribution of the 
model problem. We also found that convergence to a desired accuracy 
could be hastened by allowing the location of the dipole field sin- 
gularities at a variable distance from the reference plane. 

Quadrupole fields were also used for synthesis of solutions to the 
model problem. \^hen used in the same way as the dipole fields, we found 
that better convergence limits were achieved with about the same amount 
of computation time. In addition, we were able to develop procedures 
for the addition of one quadrupole field for each trajectory during a 
single iteration of a scheme designed to meet the constraints of all the 
trajectories, excepting the scaling trajectory, simultaneously. This 
scheme reduced by a factor of 5 (i.e., the number of trajectories con- 
sidered) , the number of iterations and the amount of computational time 
required to obtain a desired accuracy. Future possibilities for ex- 
tending this work include the development of schemes that combine the 
scaling operation and the quadrupole synthesis operation into a single 
operation which meets the constraints for all of the the trajectories. 
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This would not alter the speed of the procedures appreciably but might 
improve the scale accuracy. Also of interest are schemes which com- 
bine the dipole and quadrupole methods for potential field synthesis. 

Finally, a test potential was created in the laboratory to test 
the computer codes developed for the synthesis of unknown potentials. 
An electron-beam probe was used to collect constraining electron tra- 
jectory data. Potentials were synthesized which seemed like reason- 
able approximations considering that the biased-fin potential that 
we attempted to create in the laboratory was undoubtably altered some- 
what by the finite size of the physical apparatus. We learned that 
data selection is Important for the best determination of unknown 
potentials when applying the particular codes that we have developed. 
Therefore, the challenge in the area of dipole and quadrupole field 
synthesis lies in the determination of optimal schemes depending on 
whatever data are available. 
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APPENDIX A 


Calculating Potential Fields for 
Ideal Line Dipoles and Quadrupoles 
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To determine the potential field (j>^ existing near an ideal line 
dipole located at point (x^,y^), consider a line of point dipoles ex- 
tending to infinity in the positive and negative z-directions as in 
Figure Al. At any given distance from the x-y plane, there exists 
an electric dipole contributing a potential field. 


4)^(x,y,z^) 



cos g(x,y) 

b (x,y,Zj^) 


(A-1) 


at all points in the x-y plane (ref. 1). The electric potential at an 
arbitrary point (x,y) due to the ideal line dipole 


<t»^(x,y) 



4>^(x,y,Zj^)dZj^ 


(A-2) 


is found by integrating the potential field due to an electric dipole 
at distance from the x-y plane over all values of z^^ from - <» to + <». 
Note that 


b^(x,y,Zj^) = z^^ + (x-x^)^ + (y-y^)^ 


(A-3) 


and 


cos 8(x,y) 


[z^^ + (x-x^)^ + (y-y^)^]^^^ 


(A-4) 


Therefore, 


f '^^l 

<p,(x,y) =■ A (y-y ) 2 2 2 

d d d [ [z/ + (x-xj^ + (y-y^)'"] 


00 ^“1 


3/2' 


(A-5) 


Evaluating the integral, we find 


Ad(y-yd) 


(x-x^)^ + (y-y^)^ 


4>jj(x,y) = 


(A-6) 
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with the aid of CRC Integral Tables (ref. 8). Then, if an image line 
dipole is also present at point (x^, ~y^) » the combined electric 
potential at point (x,y) is 


(j),(x.y) = A,[- 


(y-Yd) 


iy+y,) 

+ -] . (A-7) 


.^d^ .2 . v2 ■ ^2 ^ ^ v2 

(x-x^) + (y-y^) ^y^d^ 


In a s imil ar manner, the potential (})^ existing near an ideal line 
quadrupole may be considered as an integral of the field due to a 
point quadrupole. 


A /x V 2 ■> = a sin S(x,y) cos B(x,y) 

cp ^x,y,2i; a 

^ b (x,y,z^) 


(A-8) 


located a distance from the x-y plane, over all values from - to 


<l>q(x,y) = 


(j)^(x,y,z^)dz^ 


(A-9) 


Using equations (A-3) and (A-4), along with 


sin B(x,y) = 


(x-Xd) 


222 1 / 2 ’ 
[z^ + (x-x^) + (y-t^) ] 


(A-10) 


we obtain the integral expression; 




dz. 


r 2 ^ , ^2 ^ , ^2^5/2 

-oo + (x-x^) + (y-y^) ] 


(A-11) 


Evaluating the integral, we find 


A (x-x ) (y-y ,) 


2 2 2 ’ 
[(x-x^) + (y-y^)^]^ 


(A-12) 
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where = — a^. Then, if an image line dipole is also present at point 

(x,, -y,), the combined electric quadrupole potential at point (x,y) is 
a a 


(j) (x,y) = A (x-x ) 
q q d 


(y-y,) 

[ O + o 0-J- (A-13) 

[(x-x^) + (y-y^j) ] [(x-x^) + (y+y^) 1 


For an r-0 coordinate system, consider the potential field due to 

an ideal line dipole Intersecting the r-0 plane at the point (rj,0.) 

d d 

and oriented so that the axes of all the point dipoles comprising the 
line dipole extend through points on a perpendicular line at the origin 
of the r-0 plane. This situation is shown in Figure A2. The point 
dipole located a distance from the r-0 plane contributes a potential 
field 


(|)^(r,0,Zj^) 


cos B(r,0) 
b^(r,0,zj^)’ 


(A-14) 


at all points (r ,0 ) in the r-0 plane. The electric potential due to 
m n 

the entire line dipole 


4>^(r,0) 


<j)^(r,0,Zj^)dz^ 


(A-15) 


is found by integrating equation (A-14) for all point dipoles located 
along the line (i.e. - <z^<«>). 

2 2 2 

Applying the Pythagorean Theorem, b = z^ + d , and the law of 
cosines, 

d^ = r^ + r,^ - 2rr cos(0j-0) , (A-16) 

a da 

and recognizing in Figure A2 that 
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r -r cos(0 -6) 

cos B = — r ; (A- 17) 

D 


we derive 

7 A [r -r cos(0 -0)1 

, ,2 ^ 2,3/2 ‘‘^1- 

-L (d + ) 


Then the integral may be evaluated with the result 

A [r -r cos (0-0)1 

<j).(r,0) = T~~' 2 • 

r , + r - 2rr, cos (0,-0) 

a d u 



APPENDIX B 


Fortran Routines for Dipole Synthesis of Solutions 
to The Boundary Value Problem 
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cROUTINE WIODGR.F4 

IMPLICIT REAL *8(A-K,0-Z) 

CO?-C ION /'iIVD/XB(40),YS(40),PPHI(40),IB,F,S 
CO'lC 10N/ALL/PSI(80,40) ,PHI(40) ,M,N 
COM}‘raM/MVPA/R( SO) ,T1IETA(40) 


IT=0 

CALL DATA 
:E-IAX=-1.0D38 

:giin=i.od3S 

\RLAX=-1.0D38 
Y!mi=l .0D38 
XB(IB+1)=XB(1) 

YBCIB+1)=YB(1) 

XL=0.0 

YL=0.0 

AL=0.0 

DO 120 K=1,IB 

IF ( XB ( K) . CT . >2LAX) :E1A:I=XB ( K) 

IF(XB(K) .LT.XMIM) !G!IM=XB(K) 

IF(Y3(K) .GT.Y'-LAX) Y1LA::=YB(K) 

IF(YB(K) .LT.YMIN) YMIM=YB(X) 

IF(XE(K+1) .ME.XBdO) 00 TO 100 
XL=’vL+XB ( K) *DABG (YB ( K+1 ) -YB (K) ) 

YL=YL+(Y3 (K+D+YB (X) ) *0 . 5*DABS (YB (X+1 ) -YB (K) ) 
AL=AL+D AB S ( YB ( X-!- 1 ) -YB ( K) ) 

GO TO 120 

100 IF(YB(X+1) .IIE.YB(K)) GO TO 110 

:iL =XL+ ( XB ( X+ 1 ) +XB ( X) ) *0 . 5 *D AB S ( XB ( X+ 1 ) -XB (K ) ) 
YL=YL+YB ( K) *D ABS ( XB ( K+1 ) -XB (K) ) 

AL=AL+D AB S ( XB ( K+ 1 ) -XB ( K) ) 

CO TO 120 

110 Ti\IIB= ( YB ( X+ 1) -YB ( X) ) / (XB ( K+ 1 ) -XB ( K ) ) 

XL=XL+( 1+TANB**2) **0 . 5*0 . 5* (XB (K+1) **2-X3 (K) **2) 
YL=YL+(1+1/TA?T3**2)**0.5*0.5*(YB(K+1)**2-YB(K)**2) 
AX= AL+ ( ( XB ( K+ 1 ) -XB ( K) ) ** 2+( YB ( K+ 1 ) -YB ( K) ) ** 2 ) **0 . 5 
120 CONTINUE 

X0=XL/AL 
Y0=YL/AL 
TYPE 125,X0.X0 

125 F0R1LAT(' XO=' ,D, 'Y0=' ,D) 

PHI’L\X=-1.0D3S 
PHI!-IIN=1.0D33 
PGL\X=0 . 0 
PI=4*DATM(1.D0) 

DO 290 J=1,N 
JK=0 

THETA(J)=2*PI*J/U 

160 IF(DABS(DSIN(TnETA(J))).LE.l.D-10) CO TO 170 

IF(DAB5(DC0S(TI!ETA(J))).LE.1.D-10) GO TO ISO 
TANJ=DS IN (THETA (J) ) /DCOS (TKSTA( J) ) 

GO TO 190 
JK=1 


170 
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GO TO 190 
ISO JK=-1 

190 DO 2S0 K=1,IB 

THETAl =D ATAII 2 ( YB ( K+ 1 ) -YO , XE ( K+ 1 ) -XO ) 

THETA2=D ATA1I2 ( YB ( K) -YO , XB ( K) -XO ) 

IF(T1!ETA1.GE.0.) GO TO 195 
rrfETAl=THETAl+2*PI 
195 IF(T1IETA2.GE.0.) GO TO 197 

THETA2=THETA2+2*PI 

197 IF(THETA1.GT.THETA2) GO TO 198 
IF(T11ETA(J) .LE.TUETAl) GO TO 280 
IF(THETA(J) .GT.THETA2) GO TO 280 
GO TO 888 

198 IF(THETA(J) .LT.PI) GO TO 199 
IF(THETA(J) .LE. THETAl) GO TO 280 
GO TO 888 

199 IF(THETA(J) .GT.THETA2) GO TO 280 

888 IF(DABS(XB(K)-XB(K+1)).LE.1.D-10) GO TO 230 

IF(DABS(YB(K)-YB(K+D) .LE.l.D-10) GO TO 240 

IF(JK) 210,220,200 

200 Y=Y0 

X=XB(K)-(YB(K)-Y0)*(XB(K)-XB(K+1))/(YB(K)-YB(K+D) 

GO TO 270 
210 X=X0 

Y=Y3(K)-(XB(K)-X0)*(YB(K)-Y3(K+1))/(XB(K)-XB(K+D) 

GO TO 270 

220 YN=YB ( K) * ( XB ( K) -XB ( K+ 1 ) ) / ( YB ( K) -YB ( K+1 ) ) -XB ( K) +X0+Y0/TAH J 

YD= ( XB ( K) -XB ( X+ 1 ) ) / ( YB ( K) -YB ( K+1) ) 

Y=YM/YD 

X=X0+(Y-Y0)/TAHJ 
GO TO 270 
230 X=XB(K) 

IF(JK) 280,236,232 
232 Y=Y0 

GO TO 270 

236 Y=Y0+(X-X0)*TANJ 

GO TO 270 
240 Y=YB(K) 

IF(JK) 244,246,280 
244 X=X0 

GO TO 270 

246 X=X0+(Y-Y0)/TANJ 

270 R(J)=(X**2+Y**2)**0.5 

IF ( R ( J ) . GT . RI L4X) R1L\X=R ( J ) 

PHI(J)=PPHI(K) 

IF(PHI( J) .LT.PHFilIfl) PHIMIM=PHI( J) 

IF(PHI( J) .GT.PIIF-LAX) PI1II1AX=PHI ( J) 

230 COMTIITUE 

290 CCIITINUE 

H=2*N 
IH=H+1 

DO 291 I=IN,M 
R ( I ) = ( I-il) *PwL4X/ (F-X) 
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291 CONTINUE 

DO 295 J=1,N 
DO 294 1=1, M 

PSI (I , J)=( PHIMAX-PHIMIN) *R( I) *DCOS (THETA( J) ) / ( 2*R?IAX)+(PHIMAX+PHI 
1 MIN)/ 2 

294 CONTINUE 

295 CONTINUE 
CALL ROTATE 
CO TO 301 

300 CALL SHIFT 

301 CALL SCALE 
DAN=0 . 0 

DO 310 J=1,N 

D AN=D.\N+D AG S(PHI(J)-?SI(J,J)) 

310 CONTIJIUE 

IF (DAN. LE. 0.00001) GO TO 320 
CALL TADPOLE(F.S) 

IT=IT+1 
GO TO 300 

320 PAUSE 

CALL IIIITT 

CALL DWINDO(SNGL(XNIN) .SNGLCGIAN) , SNGL ( YNIN) ,SNGL(YJLVX) ) 

CAI.L TWIMDO(0, 1023,0,780) 

CALL ERASE 

CALL VGRAPH ( XO , YO , PIIIMIN , PKIiL\X) 

PAUSE 

CALL PRINT (IT) 

STOP 

330 TYPE 340 

340 FOR2LVT(' THE POINT X0,Y0 LIES ON THE SOUNDARY') 

CALL FINITT 

STOP 

END 
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1 

2 

3 


4 

5 

6 

7 

8 

2 ? 


29 


SUBROUTINE DATA 
IMPLICIT REAL *8 (A-H.O-Z) 
COMMON/’D/D/XB(40) ,YB(40) ,PPHI (40) , IB , F , S 
COFEION/ALL /P S I( 8 0 , 4 0 ) , PHI ( 4 0 ) , M , M 
TYPE 1 


F0R1L\T(' ENTER NUMBER OF ANGULAR GRID LINES') 

ACCEPT 2,M 
FOPxMAT(I) 

TYPE 3 

FORI-tATC' ENTER NUMBER OF STRAIGHT BOUNDARY SEC'MENTS' ) 
ACCEPT 2, IB 
DO 8 K=1,IB 
TYPE 4,K 

FORMAT(' ENTER X-COORD. AT BEGINNING OF SEGMENT ',13) 
ACCEPT 5,X3(K) 

FOPJIAT(D) 

TYPE 5,K 

FORilAT(' ENTER Y-COORD. AT BEGINNING OF SEGMENT ',13) 
ACCEPT 5,Y3(K) 


TYPE 7,R 

FOR]LAT(' ENTER POTENTIAL ON BOUNDARY SEGMENT ',13) 
ACCEPT 5,PPKI(K) 

CONTINUE 


TYPE 28 

FORiL\T(' ENTER DIPOLE G PA?^\2IETER' ) 
ACCEPT 5,F 
TYPE 29 

FOPJ-LVT(' ENTER DIPOLE T PARAMETER') 

ACCEPT 5,S 

RETURN 

E!TD 
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SUBROUTINE ADPOLE(T) 

IMPLICIT REAL *8 (A-H,0-Z) 

COI-C-IOM/ALL/PSI(80,40) ,PHI(40) 

COMMON/MVPA/R(80) ,THETA(40) 

DIMENSION D(45) 

DAN=0 . 0 
DMAX=0 . 0 
DO 41 J=1,N 
D(J)=PHI(J)-PSI(J,J) 

DAN=DAM+DABS(D(J)) 

IF(DABS(D(J)) .LE.DMAX) GO TO 41 
DMAX=DABS(D(J)) 

JD=J 

41 CONTINUE 

AVDAN=DAN/N 

KL=0 

JN=JD+N 

DO 44 J=JD,JN 

IF(J.GT.N) GO TO 42 

JL=J 

GO TO 43 

42 JL=J-N 

43 IF(D(JL)/D(JD) .GT.0.0) GO TO 44 

KL=KL+1 

IF(KL.EQ.l) L=JL 
LL=JL 

44 CONTINUE 
IF(L.CT.JD) GO TO 45 
IL=JD-L 

GO TO 46 

45 IL=L-JD 

46 IF(LL.GT.JD) GO TO 805 
ILL=JD-LL 

GO TO 806 

805 ILL=LL-JD 

806 IF(IL.GT.ILL) GO TO 807 
JC=L 

GO TO 47 

807 JC=LL 

47 A=DABS(D(JD))-T*A\T)AN 

B=R( JC) *DCOS (THETA( JC) -THETA( JD) ) * ( 2*T* AVDAN-DABS (D ( JD ) ) ) - 
1 R(JD)*DABS(D(JD)) 

C=R(JD)*R(JC)*DCOS(THETA(JC)-THETA(JD))*DABS(D(JD))- 
1 R(JC)**2*T*AVDAN 

IF(B**2-4*A*C.LT.0.0) GO TO 60 
IF(B.LT.O.O) GO TO 48 
PJ)1=(2AC)/(-B-(BA*2-4*A*C)**0.5) 

RD2= (-B-( B**2-4*A*C) **0 . 5) / ( 2*A) 

GO TO 49 

48 RDl=(-3+(B**2-4*A*C)**0.5)/(2*A) 
RD2=(2*C)/(-3+(3**2-4*A*C)**0.5) 

49 IF(RD1.LE.R(JD)) GO TO 50 
RD=RD1 
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GO TO 51 

50 IF(RD2.LE.R(JD)) GO TO 61 

RD=RD2 

51 IF(RD.LT.R(JD)*(1.0+6.2S35/H)) RD=R( JD) *(6 . 2835/N+l .0) 

AD=(RD-R(JD))*D(JD) 

TYPE 515,RD,TKETA(JD) 

515 F0PJ1AT(' PJD=',D/ THETA=',D) 

M=2*IJ 

DO 53 J=1,N 
DO 52 1=1 ,M 

PSI(I,J)=PSI(I, J)+AD*(RD-R(I)*DCOS(THETA(J)-THETA(JD)))/(RD**2 
1 +R(I)**2-2*RD*R(I)*DCOS(THETA(J)-THETA( JD) )) 

52 CONTINUE 

53 CONTINUE 
RETURN 

60 TYPE 62 

62 FORIIAT(' DIPOLE RADIUS INCLUDES IMVGINARY ROOT') 

RETURN 

61 TYPE 63 

63 FOR!IAT(' DIPOLE LOCATED INSIDE BOUNDARY') 

RETURN 

END 



SUBROUTINE TADPOLE(G,T) 

IMPLICIT REAL *8 (A-R.O-Z) 

C0>E10M/ALL/PSI ( 80 , 40) , PHI (40) ,M,N 
COMMON/MVPA/R(80) ,THETA(40) 

DIMENSION D(45) 

DAN=0.0 
DMAX=0 . 0 
DO 41 J=1,N 
D(J)=PHI(J)-PSI(J,J) 

DAN=DAN+DABS(D(J)) 

IF(DABS(D(J)) .LE.DMAX) CO TO 41 
DMAX=DABS(D(J)) 

JD=J 

4 1 CONTINUE 

AVDAN=DAN/N 

KL=0 

JN=JD+N 

DO 44 J=JD,JN 

IF(J.GT.M) GO TO 42 

JL=J 

GO TO 43 

42 JL=J-N 

43 IF(D(JL)/D(JD) .GT.0.0) GO TO 44 
KL=KL+1 

IF(KL.EQ.l) L=JL 
LL«JL 

44 CONTINUE 
IF(L.GT.JD) GO TO 45 
IL=JD-L 

GO TO 46 

45 IL=L-JD 

46 IF(LL.GT.JD) GO TO 805 
ILL=JD-LL 

GO TO 806 

805 ILL=LL-JD 

806 IF(IL.CT.ILL) GO TO 807 
JC=L 

GO TO 47 

807 JC=LL 

47 A=DABS(D(JD))-T*AVDAN 

B=R(JC)*DCOS(THETA(JC)-THETA(JD) )*(2*T*AVDAN-DABS(D(JD) ))- 
1 R(JD)*DABS(D(JD)) 

C=R( JD ) *R ( JC ) *DCOS (THETA( JC ) -THETA ( JD ) ) *DABS (D ( JD ) ) - 
1 R(JC)**2*T*AVDAN 

IF(B**2-4*A*C.LT.O.O) GO TO 60 
IF(B.LT.O.O) GO TO 48 
RDl=(2*C)/(-B-(B**2-4*A*C)**0.5) 
RD2=(-B-(B**2-4*A*C)**0.5)/(2*A) 

GO TO 49 

48 RD1=(-B+(B**2-4*A*C)**0.5)/(2*A) 
RD2=(2*C)/(-B+(B**2-4*A*C)**0.5) 

49 IF(RD1.LE.R(JD)) GO TO 50 
RD=RD1 
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GO TO 51 

50 IF(RD2.LE.R(JD)) GO TO 61 

RD=RD2 

51 IF(RD.LT.R(JD)*(1.04G/M)) RD=R(JD)*(G/N+1 .0) 

AD=(RD-R(JD))*D(JD) 

TYPE 515,RD,THETA(JD) 

515 FORMAT(' RD=',D,' THETA=',D) 

M=2*N 

DO 53 J=1,N 
DO 52 1=1, M 

PSI(I,J)=PSI(I,J)+AD*(RD-R(I)*DCOS(THETA(J)-THETA(JD)))/(RD**2 
1 +R(I)**2-2*RD*R(I)*DCOS(THETA(J)-THETA(JD))) 

52 CONTINUE 

53 CONTINUE 
RETURN 

60 TYPE 62 

62 FORMAT (' DIPOLE RADIUS INCLUDES IMAGINARY ROOT') 

RETURN 

61 TYPE 63 

63 FORHAT(' DIPOLE LOCATED INSIDE BOUNDARY') 

RETURN 

END 
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SUBROUTINE SHIFT 
IMPLICIT REAL*8(A-H,0-Z) 

COMHON/ALL/PSI(80,40),PHI(40),M,N 
DIMENSION QPSK40) 

DON=0.0 
DO 85 J=1,N 

DON=DON+PHI ( J ) -P S I ( J , J ) 

85 CONTINUE 
ADON=DON/N 
DO 88 1=1 ,M 
DO 86 J=1,N 
QPSI(J)=PSI(I,J) 

86 CONTINUE 

DO 87 J=1,N 

PSI(I,J)=qpsi(J)+ADON 

87 CONTINUE 

88 CONTINUE 
RETURN 
END 

49 IF(RD1.LE.R(JD)) GO TO 50 

RD=RD1 

GO TO 51 

50 IF(PJ)2.LE.R(JD)) GO TO 61 

RD=RD2 

51 IF(RU.LT.R(JD)*(1.0+6.2835/N)) RD=R( JD) * (6 . 2835/N+l . 0) 

AD=(RD-R(JD))*D(JD) 

TYPE 515,RD,THETA(JD) 

515 FOR.’tAT(' RD=',D/ THETA=',D) 

M=2*N 

DO 53 J=1,N 
DO 52 1=1, II 

PSI(I,J)=PSI(I, J)+AD*(RD-R(I)*DCOS(THETA(J)-THETA(JD)))/(RD**2 
1 +R(I)A*2-2 ARD*R ( I ) *DCO S ( THETA ( J ) -THETA ( JD ) ) ) 

52 CONTINUE 

53 CONTINUE 
RETURN 

60 TYPE 62 

62 FOP2IAT(' DIPOLE RADIUS INCLUDES II-AGINARY ROOT') 

RETURN 

61 TYPE 63 

63 FORILATC' DIPOLE LOCATED INSIDE BOUNDARY') 

RETURN 

ElTD 

SUBROUTINE VGPu\PH(XO,YO ,PHIHIM, PHIMAIO 
IMPLICIT REAL *8 (A-H,0-Z) 

COMMON/M\TD/X3(40) ,YB(40) ,PPHI(40) ,IB,S 
COMMON/ALL/PSI(80,40) ,PHI(40) ,M,N 
COMMON/HVPA/R(30) ,THETA(40) 

DIMENSION X(190) ,Y(190) 

C.U.L MOVEA( SNCL (XB ( IB) ) , SNGL(YB ( 13) ) ) 

DO 500 K=1,IB 

C.\LL DRAUA( SNGL ( >EB ( K) ) , SNCL ( YB ( IC) ) ) 

500 CONTINUE 



P1=4*DATAN(1.D0) 

MR=N+1 

Mi-I=N+2 

DO 595 K=l,39 

V=PHIMIM+K* (PHII'IAX-PIIIIIIN) /40 
L=0 

DO 515 J=1,N 

DO 510 I=M1I,H 

IF(R(I) .GT.R(J)) CO TO 515 

IF((PSI(I-1,J)-V)*(PSI(I,J)-V) .GT.0.0) GO TO 510 
L=L+1 

RR=R(I)+(V-PSI(I, J))*(R(I)-R(I-1))/(PSI(I,J)-PSI(I-1,J)) 
X(L) =RR*DCOS (T11ETA( J) ) 

Y ( L) =RR*D S IN ( THET A ( J ) ) 

510 ' CONTINUE 

515 CONTINUE 

THETA(N+1 ) =THETA( 1 ) 

DO 525 I=MR,N 
PSI(I,N+1)=PSI(I, 1) 

DO 520 J=2,MR 
JJ=J 

IF(JJ.EQ.N+1) JJ=1 

IF(R(I) .GT.R(JJ)) GO TO 520 

IF(R(I) .GT.R(J-l)) GO TO 520 

IF((PSI(I,J-1)-V)*(PSI(I,J)-V).GT.0.0) GO TO 520 
L=L+1 

TT=THETA(J)+(V-PSI(I,J))*2*PI/((PSI(I,J)-PSI(I,J-1))*N) 

X(L)=R(I)*DCOS(TT) 

Y(L)=R(I)*DSIH(TT) 

520 CONTINUE 

525 CONTINUE 

LI-IAX=L 
DG=0.0 

DO 570 L=l,LILAIv 

IF((X0-X(L))**2+(Y0-Y(L))**2.LT.DG) GO TO 570 

DG=(X0-X(L) ) **2+(Y0-Y(L) ) **2 

LY=L 

570 CONTINUE 

:ci=x(i) 

YY=Y(1) 

X(1)=X(LY) 

Y(1)=Y(LY) I 

X(LY)=XX 

Y(LY)=YY 

CALL MOVEA(SNGL(X(i)),SNGL(Y(l))) 

LL=1 

535 LL=LL+1 

DSQ=1.0D38 
DO 540 L=LL,LMA.X 

DSQL=(X(L)-X(LL-1))**2+(Y(L)-Y(LL-1))**2 

IF(DSQL.GT.DSQ) GO TO 540 

DSQ=DSQL 

LDMIM=L 
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540 CONTINUE 

XX=X(LL) 

YY=Y(LL) 

X(LL)=X(LDMIN) 

Y(LL)=Y(LDMIN) 

X(LDMIM)=XX 

Y(LDM1N)=YY 

CALL DRAWA(SNGL(X(LL) ) ,SNGL(Y(LL) ) ) 
IF(LL.EQ.LMAX) GO TO 595 
GO TO 535 
595 CONTINUE 

RETURN 
END 
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SUBROUTINE SCALE 
IMPLICIT REAL *8 (A-H.O-Z) 

COMI-10N/ALL/PSI(80,40) ,PHI(40) ,M,N 

AVPHI=0.0 

AVPSI=0.0 

DO 80 J=1,M 

AVPHI=AVPHI+PHI ( J ) 

AVPSI=AVPSI+PSI(J,J) 

80 CONTINUE 
AVPHI=AVPHI/N 
AVPSI=AVPSI/N 
ABPHI=0.0 
ABPS 1=0.0 

DO 81 J=1,M 

ABPHI=ABPHI+DA3S(PHI(J)-AVPHI) 

ABPSI=ABPSI+DABS(PSI(J,J)-AVPSI) 

8 1 CONTINUE 
M=2*N 

DO 83 J=1,N 
DO 82 1=1 ,M 

IF(ABPSI.EQ.n.O) CO TO 82 

PSI(I, J) = (PSI(I,J)_,'^VPSX)*A3PHI/ABPSI+AVPSI 
■ 82 CONTINUE 

83 CONTINUE 

RETURN 
END 

49 IF(RD1.LE.R(JD)) GO TO 50 

RD=RD1 

GO TO 51 

50 IF(RD2.LE.R(JD)) CO TO 61 

PJ)=PJ)2 

51 IF(RD.LT.R(JD)*(1.0+6.2835/N)) RD=R( JD)*(6.2835/N+1 .0) 

AD=(PJD-R(JD))*D(JD) 

TYPE 515,RD,TT!ETA(JD) 

515 FORIU'.TC' ?J5=',D,' THETA=',D) 

M=2*N 

DO 53 J=1,N 
DO 52 1=1, M 

PSI(I,J)=PSI(I,J) +j\D* (PJD-R ( I ) *DCOS ( THETA ( J ) -THETA ( JD ) ) ) / (RD* * 2 
1 +R(I)**2-2*RD*R(I)*nCOS(THETA(J)-THETA(JD) ) ) 

5 2 CONTINUE 

53 CONTIIHJE 

RETURN 

60 TYPE 62 

62 FORI‘L\T(' DIPOLE RADIUS INCLUDES IMAGINARY ROOT") 

RETURN 

6 1 TYPE 6 3 

63 F0P1U\T(' DIPOLE LOCATED INSIDE BOUNDARY') 

RETURN 

END 

SUBROUTINE VGRAPI! ( NO , YO , PHIMIN , PlIINAlO 
IMPLICIT REAL *3 (A-H,0-Z) 

CO:G!ON/:!\U)/XB(40) , YB (40) , PPIII (40) , IB , S 
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C0>C-IOn/ALL/PSI(8O,40) , PHI (40) ,H,H 
COMMON /MVPA /R (80), THETA (40) 

DIMENSION X(190) ,Y(190) 

CALL MO VE A ( SNGL ( XB ( IS ) ) , SNGL ( YB ( IB ) ) ) 

DO 500 K=1,IB 

CALL DRAMA ( SNGL (XB(K) ) ,SNGL(YE(K) ) ) 

500 CONTINUE 

PI=4*DATAN(1.D0) 

MR=N+1 

>[N=N+2 

DO 595 K=l,39 

V=PHIMIK+H* (PHDIAX-PHIHIN) /40 
L=0 

DO 515 J=1,N 
DO 510 I=IC1,M 
IF(R(I).GT.R(J)) GO TO 515 

IF((PSI(I-1,J)-V)*(PSI(I,J)-V).GT.0.0) GO TO 510 
L=L+1 

RR=R(I)+(V-PSI(I,J))*(R(I)-R(I-1))/(PSI(I,J)-PSI(I-1,J)) 
X(L)=RR*DCOS(THETA( J)) 

Y ( L ) =RR*D S IN ( THETA ( J ) ) 

510 CONTINUE 

515 CONTINUE 

THETA(N+1)=THETA( 1) 

DO 525 I=MR,M 
PSI(I,N+1)=PSI(I, 1) 

DO 530 J=2,MR 
JJ=J 

IF(JJ.EQ.N+1) JJ=1 

IF(R(I) .GT.R(JJ)) GO TO 520 

IF(R(I) .GT.R(J-l)) GO TO 520 

IF( (PSI(I, J-1)-V)*(PSI(I,J)-V) .GT.0.0) GO. TO 520 
L=L+1 

TT=THETA(J)+(V-PSI(I,J))*2*PI/((PSI(I,J)-PSI(I, J-1))*N) 
X(L)=R(I)*DCOS(TT) 

Y(L)=R(I)*DSIN(TT) 

520 CONTINUE 

525 CONTINUE 

LMAX=L 
DG=0.0 

DO 570 l=i,l;iax 

IF((X0-X(L))**2+(Y0-Y(L))**2.LT.DG) CO TO 570 

DG=(X0-X(L))**2+(Y0-Y(L))**2 

LY=L 

570 CONTINUE 

XX=X(1) 

YY=Y ( 1 ) 

X(1)=X(LY) 

Y(1)=Y(LY) 

X(LY)=XX 

Y(LY)=YY 

C.\LL MOVEA( SNGL ( X ( 1 ) ) , SNGL ( Y ( 1 ) ) ) 

LL=1 
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535 LL=LL+1 

DSQ=1.0D33 
DO 540 L=LL,LHAX 

DSQL=(X(L)-X(LL-1))**2+(Y(L)-Y(LL-1))**2 

IF(DSQL.GT.DSQ) GO TO 540 

DSQ=DSQL 

LDMIN=L 

540 COtlTIMUE 

XX=X(LL) 

YY=Y(LL) 

X(LL)=X(LDMIN) 

Y(LL)=Y(LDMIM) 

X(LDMIN)=XX 

Y(LDIIIN)=YY 

CALL DRAUA(SNGL(X(LL) ) ,SnGL(Y(LL) ) ) 
IF(LL.EO.LMAX) CO TO 595 
GO TO 535 
595 CONTINUE 

RETURN 
END 



APPENDIX C 


Fortran Routines for Generating the Biased Fin 
Potential and the Associated Electron Trajectory Data 
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cROUTINE TEST.F4 

IMPLICIT REAL *8 (A-II.O-Z) 

COMPLEX Z,ZZ 

DIMENSION PHI(82,22) ,XX(200) ,YY(200) ,X(82) ,Y(22) 


PI=4*DATAN(1.0D0) 

TYPE 10 

10 FOPC>IAT(' ENTER X DIMENSION=' , $) 

ACCEPT 20,DX 
20 FORl’IAT(D) 

TYPE 30 

30 FORMAT (' ENTER Y DIMENSION=' , $) 

ACCEPT 20.DY 
X(l)=-DX/30 
DO 40 1=2,82 
X(I)=(I-2)*DX/30 
40 CONTINUE 

Y(l)=-DY/20 
DO 50 J=2,22 
Y(J)=(J-2)*DY/20 
50 CONTINUE 

S=PI/1.8 

PSI=DLOG(-DEXP(S)+DSQRT(l+DEXP(2*S))) 

TYPE 55,PSI . 

55 FORIIAT(D) 

DO 70 1=22,41 
DO 60 J=2,22 
S=PI*(1-Y(J)/DY)/1.8 
T=pi*(2*X(I)/DX-0.5) 

Z=CMPLX ( SNGL ( S ) , SNGL (T) ) 

ZZ=CSQRT ( 1-CEXP ( 2*Z ) ) 

D=REAL(ZZ) 

E=AIMAG(ZZ) 

ZZ=CEXP(Z) 

DD=-AIMAG(ZZ)‘ 

EE=REAL(ZZ) 

ZZ=CMPLX(SMGL(D+DD) ,SNGL(E+EE) ) 

ZZ=CLOG(ZZ) 

PHI(I, J)=-REAL(ZZ)-PSI 
60 CONTINUE 

70 CONTINUE 

DO 75 J=2,22 
S=PI*(1-Y(J)/DY)/1.3 

PHI(42,J)=DLOG(-DEXP(S)+DSORT{l+DEXP(2*S)))-PSI 
75 . CONTINUE 

DO 90 K=2,21 
DO 80 J=2,22 
PHI(K,J)=PHI(44-K, J) 

80 CONTINUE 

90 CONTINUE 

DO 110 L=3,42 
DO 100 J=2,22 
PHK40+L, J)=PHI(L, J) 



112 


100 CONTINUE 

110 CONTINUE 

CALL INITT 

CALL DWINDO(SNGL(X(2) ) ,SNGL(X(82) ) ,SNGL(Y(2) ) ,SNGL(Y(22) ) ) 
II=INT(SNGL(780*DY/DX)) 

CALL TWINDO(0, 780,0,11) 

CALL MO VE A ( SNGL ( X ( 2 ) ) , SNGL ( Y ( 2 ) ) ) 

CALL DRAWA(SNGL(X(82)) ,SNGL(Y(2))) 

PHIMIN=1.0D38 
PHIMAX=-1.0D38 
DO 140 1=2,82 
DO 130 J=2,22 

IF(PHI(I,J).GT.PHIMIN) GO TO 125 
PHIMIN=PHI(I,J) 

125 IF(PHI(I,J) .LT.PHIMAX) GO TO 130 
PHIMAX=PHI(I,J) 

130 CONTINUE 

140 CONTINUE 

TYPE 290,PHIMIN,PHIMAX 

290 FORMAT(' PHIMIN=' ,D12.6,' PHIMAX=' ,D12.6) 

K=0 

142 K=K+1 

V=PHIMIN+K*0 . 2 
L=0 

DO 160 J=2,22 
DO 150 1=3,82 

IF((PHI(I-1,J)-V)*(PHI(I,J)-V).GE.0.0) GO TO 150 
L=L+1 

XX(L)=X(I-1)+(X(I)-X(I-1))*(V-PHI(I-1,J))/(PHI(I,J)-PHI(I-1,J)) 

YY(L)=Y(J) 

150 CONTINUE 

160 CONTINUE 

DO 180 1=2,82 
DO 170 J=3,22 

IF((PHI(I,J-1)-V)*(PHI(I,J)-V).GE.0.0) GO TO 170 
L=L+1 

XX(L)=X(I) 

YY(L)=Y(J-1)+(Y(J)-Y(J-1))*(V-PHI(I,J-1))/(PHI(I,J)-PHI(I,J-1)) 
170 CONTINUE 

180 CONTINUE 

LM=L 

XMIH=X(82) 

DO 190 L=1,LM 

IF(XX(L).GT.XMIN) GO TO 190 
XMIN=XX(L) 

LMIN=L 

190 CONTINUE 

XS=XX(1) 

YS=YY(1) 

XX(1)=XX(LMIN) 

YY(1)=YY(LMIN) 

XX(LMIN)=XS 

YY(LMIN)=YS 
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CALL MO VEA( SNGL ( XX ( 1 ) ) , SNGL ( YY ( 1 ) ) ) 

LL=1 

195 LL=LL+1 

DSQ=1.0D38 
DO 210 L=LL,LM • 

DSQL=(XX(L)-XX(LL-1) )**2+(YY(L)-YY(LL-l) )**2 

IF(DSQL.GT.DSQ) GO TO 210 

DSQ=DSQL 

LDM=L 

210 CONTINUE 

XS=XX(LL) 

YS=YY(LL) 

XX(LL)=XX(LDM) 

YY(LL)=YY(LDM) 

XX(LDM)=XS 

YY(LDM)=YS 

CALL DRAWA( SNGL ( XX(LL) ), SNGL (YY(LL) ) ) 
IF(LL.EQ.LM) GO TO 220 
GO TO 195 

220 IF(V+0.2.LE.PH1MAX) GO TO 142 

PAUSE 

OPEN (UNIT=1 , FILE= ' POTL.DAT' ) 

IWITE (1,20) X(l) 

WRITE (1,20) Y(l) 

DO 240 1=2,82 
VmiTE(l,20) X(I) 

DO 230 J=2,22 
IF(I.GT.2) CO TO 222 
WRITE! 1,20) Y(J) 

222 WRITE(1,20) PH1(I,J) 

225 FORMAT (D 12. 8) 

230 CONTINUE 

240 CONTINUE 

STOP 
END 
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cROUTINE DATA.F4 

IMPLICIT REAL *8(A-H,0-Z) 

COMMON/COORD /PHI (8 2, 22) ,XX(82) ,YY(22) ,DX,DY 
COMMOH/DATA/XO(50) ,Y0(5O) ,VX0(50) ,VY0(50) ,XT(50) ,YT(50) 
COMMON/SAVE/X(250) ,Y(250) ,FX(250) ,FY(250) ,DT(250) ,KT 
OPEN (UNIT=1 ,FILE='POTL. DAT') 

READ (1,50) XX(1) 

READ(1,50) YY(1) 

DO 200 1=2,82 
READ (1,50) XX(I) 

50 FORMAT (D) 

DO 100 J=2,22 
IF(I.GT.2) GO TO 80 
READ(1,50) YY(J) 

80 READ(1,50) PHI(I,J) 

100 CONTINUE 

200 CONTINUE 

CALL INITT 

CALL DWItroO(SNGL(XX(l)),SNGL(XX(82)) ,SNGL(YY( 1) ) ,SNGL(YY(22) ) ) 
DX=XX(82)-XX(2) 

DY=YY(22)-YY(2) 

II=INT ( SNGL ( 7 80*DY/DX) ) 

CALL TWINDO(0, 780,0,11) 

K=0 

300 K=K+1 

X0(K)=XX(13+K*9) 

Y0(K)=0.0 

CALL MOVEA( SNGL (XO (K) ) , SNGL ( YO (K) ) ) 

VX0(K)=1.0 

VY0(K)=7.35D5 

KT=0 

DT(1)=DY/(100*VY0(K)) 

400 KT=KT+1 

CALL DVOGEL(K) 

CALL DRAUA ( S NGL ( X ( KT+1 ) ) , SNGL ( Y ( KT+ 1 ) ) ) 

IF(Y(KT+1) .LE.0.0) GO TO 500 
GO TO 400 

500 XT(K)=X(KT+1) 

YT(K)=0.0 

IF(K.EQ.5) GO TO 600 
GO TO 300 
600 K=K+1 

X0(K)=XX(21) 

Y0(K)=0.0 

CALL MOVEA( SNGL (X0(K) ) , SNGL(Y0 (K) ) ) 

VX0(K)=5.5D5 

VY0(K)=5.0D5 

KT=0 

DT(1)=DY/(100*VY0(K) ) 

700 KT=KT+1 

CALL DVOGEL(X) 

CALL DRAUA( SNGL (X(KT+1 ) ) , SNGL ( Y (KT+1 ) ) ) 

IF(Y(KT+1) .LE.0.0) GO TO 800 
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GO TO 700 

800 XT(K)=X(KT+1) 

YT(K)=0.0 

OPEN (UNIT=1 ,FILE=' TRAJ.DAT' ) 

IL=6 

l^fRITE( 1,1500) IL 

DO 1400 K=l,6 

TOITE( 1,1600) X0(K),VX0(K) 

17RITE( 1,1600) Y0(K),VY0(K) 

VJRITE( 1,1600) XT(K),YT(K) 

1400 CONTINUE 

1500 FORIIAT(I) 

1600 F0R1IAT(2D) 

CALL VPLOT 

STOP 

END 

SUBROUTINE NEinON ( X, Y , EX, EY , PSI , IK) 

IMPLICIT REi\L *8(A-H,0-Z) 

COMMON/COORD/PHI (82, 22) ,:iX(82) ,YY(22) ,DX,DY 
DIMENSION FM(4,4,4) ,FMN(4,4,4) ,XM(4) ,1C'I1(4) ,YN(4) ,YN1(4) 
IF(X.GT.XX(82)) GO TO 220 
IF(X.LT.XXd)) GO TO 220 
IF(Y.GT.YY(22)) CO TO 220 
IF(Y.LT.YYd)) GO TO 220 
HX=DX/30.0 
HY=DY/20.0 
1=0 

10 1 = 1+1 

IF(X.GT.XXd)) GO TO 10 
IF(I.EQ.2) GO TO 12 
IB=I-2 
GO TO 14 
12 IB=I-1 

14 IF(I.EQ.82) GO TO 16 

IE=I+1 
GO TO 18 
16 IE=I 

18 J=0 

20 J=J+1 

IF(Y.GT.YY(J)) GO TO 20 
IF(J.EQ.2) GO TO 25 
J3=J-2 
GO TO 30 
25 JB=J-1 

30 IF(J.EQ.22) GO TO 35 

JE=J+1 
CO TO 40 
35 JE=J 

40 MM=IE-IB+1 

:'2‘I(1) = 1 .0 

x:iid)=o.o 

DO 70 M=2,HM 
IM=IC+M-2 
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XM ( M ) =X!-I ( M- 1 ) * ( X-XX (III)) 

XM1(M)=0.0 
DO 60 1=2 ,!I 
Xl=1.0 
DO 50 K=2,M 
IF(K.EQ.I) GO TO 50 
X 1 =X 1 * ( X-XX ( IB+K- 2 ) ) 

50 CONTI>IUE 

XM1(H)=XI-11(M)+X1 
60 COMTIIIUE 

70 CONTINUE 

NN=JE-JB+1 
YN(1)=1.0 
YN1(1)=0.0 
DO 100 N=2,NN 
JN=JB+N-2 

YN (N ) =YN (N- 1 ) * ( Y-YY ( JM ) ) 

YN1(N)=0.0 

DO 90 J=2,M 

Yl=1.0 

DO 80 L=2,M 

IF(L.EQ.J) GO TO 80 

Yl=Yl*(Y-YY(JB+L-2)) 

80 CONTINUE 

YN1(N)=Y1U(N)+Y1 
90 CONTINUE 

100 CONTINUE 

DO 110 M=l,!Dt 
DO 105 N=1,NN 

FM (1 ,H , N ) =PHI ( IB+M- 1 , JS+N- 1 ) 

105 CONTINUE 

110 CONTINUE 

DO 140 L=2,NN 
DO 130 N=L,NN 
DO 120 M=1,MM 

F!ia,M,N) = (FM(L-l,M,M-l)-FMa-l,M,M))/(HY*(l-L)) 

120 CONTINUE 

130 CONTINUE 

140 CONTINUE 

DO 160 M=1,ICI 
DO 150 N=1,NN 
FIIN(1,M,N)=FM(N,M,N) 

150 CONTINUE 

160 CONTINUE 

DO 190 K=2,^u•I 
DO 180 M=K,b21 
DO 170 N=1,NN 

FMN ( K , M , N ) = ( FMN ( K- 1 , M- 1 , N ) -FHN ( K- 1 , M , N ) ) / ( HX* ( 1 -K ) ) 
170 CONTINUE 

ISO CONTINUE 

190 CONTINUE 

EX=0.0 
EY=0.0 
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PSI=0.0 
DO 210 M=1,NN 

DO 200 M=1,1!M 

EX=EX-XM1 (M) *YK (N) *F1!N(H,M,M) 

EY=EY-YN1 (N) *XM(M) *FMM(M,M,N) 
PSI=PSI+XM(M)*YN(N)*F1!N(M,M,N) 

200 CONTINUE 

210 CONTINUE 

GO TO 230 
220 IK=-1 

230 RETURN 

END 

SUBROUTINE DVOGEL(K) 

IMPLICIT REAL *8(A-I1,0-Z) 

COMMON/COORD/PIII(82,22),XX(82) ,YY(22) ,DX,DY 
COMMON/DATA/XO(50) ,Y0(5O) ,VX0(50) ,VY0(50) ,XT(50) ,YT(50) 
COMMON/SAVE/X(250) ,Y(250) ,FX(250) ,FY(250) ,DT(250) ,KT 
0M=1.76D11 

DIM=DSQRT(DX**2+DY**2) 

IK=0 

IF(KT.GT.l) GO TO 200 
VX=VX0(K) 

VY=VY0(K) 

X(1)=X0(K) 

Y(1)=Y0(K) 

100 CALL NEUTON(X(KT),Y(KT),EX,EY,PSI,IK) 

IF(IK.NE.O) GO TO 400 

AX=QM*EX 

AY=QM*EY 

XN=X (KT ) - ( 0 . 5 * VX-0 . 1 2 5*AX*DT ( KT ) ) *DT ( KT ) 

YH=Y ( KT ) - ( 0 . 5 * VY-0 . 1 2 5 *AY*DT ( KT ) ) *DT ( KT ) 

C.\LL NEWTON ( XH , YH , EX , EY , PS I , IK) 

IF(IK.NE.O) GO TO 400 

a:<h=qm*ex 

AYH=QM*EY 

200 XII=X(KT) + (0.5*VX+(4.0*AX-A:ai)*DT(KT)/24.0)*DT(KT) 

YH=Y ( KT ) + ( 0 . 5 * VY+ ( 4 . 0*AY-AYH ) *DT ( KT ) / 2 4 . 0 ) *DT ( KT ) 

CALL NEWTON ( XH , YF. , EX , EY , PS I , IK) 

IF(IK.NE.O) GO TO 400 

AIffl=QM*EX 

AYH=QM*EY 

X(KT+1)=X(KT)+(VX+(AX+2.0*AXH)*DT(KT)/6.0)*DT(KT) 

Y ( KT+1 ) =Y ( KT ) + ( VY+( AY+2 . 0*AY!i) *DT ( KT ) / 6 . 0 ) *DT ( KT ) 

AAX=AX 

AAY=AY 

CALL NEWTON(X(KT+l) ,Y(KT+1) ,EX,EY, PSI , IK) 

IF(IK.NE.O) GO TO 400 

AX=QM*EX 

AY=QM*EY 

FX(KT)=EX 

FY(KT)=EY 

A 1 S Q= Ai\X* * 2+ AAY * * 2 

DELD =D S QRT ( ( X ( KT+ 1 ) -X { KT ) ) * * 2+ ( Y ( KT+ 1 ) - Y ( KT ) ) * * 2 ) 
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IF(DELD.GT. DIM/80) GO TO 320 
IF(DELD.LT. DIM/400) GO TO 330 
IF(A1SQ.LT. l.D-6) GO TO 300 
280 VX=VX+DT(KT)*(AX+4.0*AXH+AAX)/6.0 

VY=VY+DT ( KT ) * ( AY+4 . 0*AYH+AAY ) / 6 . 0 
DT(KT+1)=DT(KT) 

GO TO 440 

300 V1SQ=VX**2+VY**2 

DVSQ=A1SQ*DT(KT) **2 
IF(V1SQ.GT.DVSQ*1.D4) GO TO 280 
IK=1 

GO TO 410 

320 DT(KT)=DT(KT)/1.5 

GO TO 100 

330 DT(KT)=DT(KT)*1.5 

GO TO 100 
400 TYPE 420 

GO TO 440 
410 TYPE 430 

420 FORMAT (' ELECTRON HAS ESCAPED POTENTIAL REGION 

1 INCLUDED IN X-Y COORDINATE GRID') 

430 FORM/\T(' PARTICLE VELOCITY AND iVCCELEPATION BOTH 

1 EQUAL ZERO') 

440 RETURN 

END 

SUBROUTINE VPLOT 
IMPLICIT REAL *3 (A-H.O-Z) 

C0!G10N/C00PJD/PHI(S2,22) ,XX(82) ,YY(22) ,DX,DY 
DIMENSION X(180) ,Y(180) 

CALL INITT 

CALL DUINDO(SNGL(XX(D) ,SNCLao:(S2)) ,SNGL(YY(D) ,SNGL(YY(22))) 
II=INT(SNGL(7S0*DY/DX)) 

CALL TUINDO(0, 780, 0,11) 

CALL MOVEA( SNGL CLX ( 2 ) ) , SNGL (YY ( 2 ) ) ) 

CALL DPJVUA ( SNGL ( XX ( 8 2 ) ) , SNGL ( YY ( 2 ) ) ) 

PHIMIN=1 .0D38 
PHHIAX=-1.0D38 
DO 10 1=2,82 
DO 5 J=2,22 

IF(PHI(I,J) .CT.PHIMIN) GO TO 2 
PHIMIN=PHI(I, J) 

2 IF(PHI(I,J) .LT.PHIMAX) GO TO 5 
PHIMAX=PHI(I,J) 

5 CONTINUE 

10 CONTINUE 

K=0 

12 K=K+1 

V=PHIMIN+K*0 . 2 
L=0 

DO 20 J=2,22 
DO 15 1=3,82 

IF((PHI(I-1,J)-V)*(PHI(I,J)-V).GE.0.0) GO TO 15 
L=L+1 
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X(L)=XX(I-1)+(XX(I)-XX(I-1))*(V-PHI(I-1,J))/(PHI(I, J)- 
1 PHI(I-1,J)) 

Y(L)=YY(J) 

15 CONTIMUE 
20 CONTIllUE 

DO 30 1=2,82 
DO 25 J=3,22 

IF((PHI(I, J-1)-V)*(PHI(I,J)-V) .CE.0.0) CO TO 25 

L=L+1 

X(L)=XX(I) 

Y(L)=YY(J-1)+(YY(J)-YY(J-1))*(V-PHI(I,J-1))/(PHI(I, J)- 
1 PHI(I,J-D) 

25 CONTINUE 
30 CONTINUE 
LMAX=L 
XMIN=XX(82) 

DO 35L=1,LMAX 

IF(X(L) .GT.XMIN) CO TO 35 

XMIN=X(L) 

U[IN=L 

35 CONTINUE 

XS=X(1) 

YS=Y(1) 

X(1)=X(LMIN) 

Y(1)=Y(LIIIN) 

X(LIIIM)=XS 

Y(LMIN)=YS 

CALL HOVEA( SMGL ( X ( 1 ) ) , SNGL ( Y ( 1 ) ) ) 

LL=1 

40 LL=LL+1 

DSQ=1.0D3S 
DO 45 L=LL,LMAX 

DSQL=(X(L)-X(LL-1))**2+(Y(L)-Y(LL-1))**2 

IF(DSQL.GT.DSQ) GO TO 45 

DSQ=DSQL 

LDMIN=L 

45 CONTINUE 

XS=X(LL) 

YS=Y(LL) 

X(LL)=X(LDMIN) 

Y(LL)=Y(LDMIN) 

X(LDMIM)=XS 

Y(LDMIN)=YS 

CALL DPJll JA ( SNGL ( X (LL ) ) , SNGL ( Y ( LL ) ) ) 

IF(LL.EQ.LMAX) CO TO 50 
GO TO 40 

50 IF(V+0.2.LE.PHIMAX) GO TO 12 

CALL FINITT(800,500) 

RETURN 

END 
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APPENDIX D 

Fortran Routines for Dipole Synthesis of Solutions 
to The Model Problem 
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cROUTINE PTE11T.F4 

IMPLICIT REAL *8(A-H,0-Z) 

CO!E!ON/COORD/PHI(82,32) ,XX(32) ,YY(32) ,DX,DY 
CO>E'fON/DATA/X0(50),Y0(50),VX0(50),lT0(50) ,XT(50) , YT(50) 
CO?-EION/SAVE/X(110) ,Y(110) ,DT(250) ,KT 

DIMENSION XS(110),YS(110),SX(110),SY(110) ,TD(110) ,DDX(50) 

1 ,YN(50) ,XN(50) 

EXTERNAL TEKHAN.PLTHAN.SCPHAN 
CALL INITTG 
TYPE 1 

1 FORI-IAT(' ENTER X AflD Y DIMENSIONS OF POTENTIAL REGION', $) 
ACCEPT 2,DX,DY 

2 FOPvMAT(2D) 

TYPE 19 

19 FOR!LAT(' ENTER 1 TO INPUT DATA FROM FILE TRAJ.DAT, ENTER 0 
I TO INPUT DATA MAITOALLY' , $) 

ACCEPT 4, NT 
IF(NT.EQ.l) GO TO 140 
TYPE 3 

3 FORMAT(' ENTER NUMBER OF ELECTRON TPJ^JECTORIES' , $) 

ACCEPT 4,KI<; 

4 FORMAT(I) 

DO 8 K=l,iai 
TYPE 5 

5 FORMAT (' ENTER INITIAL X COORDINATE AND VELOCITY',$) 

ACCEPT 2,X0(K) ,VX0(K) 

TYPE 6 

6 FORMAT (' ENTER INITIAL Y COORDINATE AND VELOCITY ',$) 

ACCEPT 2,Y0(K) ,VYO(K) 

TYPE 7 

7 FOPvMATC' ENTER FINAL X AND Y COORDINATES' , $) 

ACCEPT 2,XT(K) ,YT(K) 

type 67,k,xO(k) ,xt (k) 

67 format(12,2dl2.6) 

8 CONTINUE 

9 EQ=1.60210D-19 
EM=9.1091D-31 
QM=EQ/EM 

A0= 2 *VXO ( KK) * VYO ( KK) / ( QM* ( XT ( KK) -XO ( KK) ) ) 

TYPE 15, AO 

15 FORMAT (' AO=',D) 

XX(l)=-DX/80 
XX(2)=0.0 
DO 20 1=3,82 
XX(I)=XX(I-l)+DX/80 

20 CONTINUE 
YY(l)=-DY/20 
YY(2)=0.0 

DO 30 J=3,32 
YY(J)=YY(J-l)+DY/20 
30 CONTINUE 

32 DO 50 1=1,82 

DO 40 J=l,32 
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PHI(I, J)=AO*YY(J) 

40 CONTINUE 

50 CONTINUE 

PAUSE 
IT=0 

60 DD=0.0 

D=0.0 

TYPE 65, IT 

65 F0RMAT(' ITERATION=' ,14) 

IT=IT+1 

DO 90 K=1,KK-1 

Y^I(K)=0.0 

KT=0 

TYPE 82,X0(K) ,Y0(K) 
DT(1)=DY/(100*'/Y0(K)) 

70 XT=KT+1 

CALL DVOGEL(K) 

IF(Y(KT+1) .LT.YK(K)) GO TO 75 
YN(K)=Y(KT+1) 

XN(K)=X(KT+1) 

75 IF(Y(KT+1) .GT.YO(K)) GO TO 70 

80 TYPE 82,X(KT+1) ,Y(KT+1) 

DDX(K)=(XT(K)-X(KT+1))/(XT(K)-X0(K)) 
82 F0R11AT(' X=',D12.4,' Y=',D12.4) 

DD=DD+DA3S ( 7H: ( K) -X ( KT+1 ) ) 
IF(D.GT.DA3S(DDX(K))) GO TO 90 
D=DABS(DDX(K)) 

XM=X(KT+1) 

roi=K 

DO 85 KI=1,KT+1 
XS(KI)=X(KI) 

YS(KI)=Y(KI) 

SX(KI)=FX(KI) 

SY(KI)=FY(KI) 

TD(KI)=DT(KI) 

85 CONTINUE 

90 CONTINUE 

TYPE 95, D 
TYPE 95,DD 
95 FORIIAT(D) 

IF(DD.LT.0.01*DX) GO TO 130 
K=KM 

180 K=K+1 

IF(K.EQ.KK) GO TO 190 
IF(DDX(K)*DDX(K-1) .GT.0.0) GO TO 180 
LP=K-KM 
GO TO 200 
190 LP=KK-KM 

200 K=KM 

210 K=K-1 

IF(K.EQ.O) GO TO 220 
IF(DDX(K)*DDX(K+1) .GT.0.0) GO TO 210 
LN=K-KH 
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GO TO 230 
220 LN=-KM 

230 IF(LP.LT.-L>i) GO TO 240 
IF(LP.GT.-LII) GO TO 235 

IF(DABS(DDX(I:M+LP)).GT.DABS(DDX(KM+LM))) go to 240 
235 L=LN 

GO TO 250 
240 L=LP 

250 IFCL.EO.-iOI) GO TO 260 

IF(L.EQ.KK-KM) GO TO 270 
XL=0 . 5 * ( XN ( K-I+L ) +Xl;I ( KL I+L- 1 ) ) 

GO TO 300 

260 XL=1.5*X1'I(1)-0.5*X1I(2) 

GO TO 300 

270 XL = 1 . 5 * XN ( KK- 1 ) - 0 . 5 *XM ( FJC- 2 ) 

300 YL=1.4*DY 

VXI 1=0.0 
DVXI 1=0.0 
VXXI=0.0 
DO 98 1=1, KT 
DVXI =D VXI I 

DVXII=2*(XS(I)-XL)*((YL-YS(I))/((XS(I)-XL)**2 

1 +(YL-YS(I))**2)**2-(YS(I)+YL)/((XS(I)-XL)**2 

2 +(YS(I)+YL)**2)**2) 

VXI = VXI I 

\GCI I=VXI+0 . 5* (D^/XI I+DVXI ) *TD ( I ) 

vxxi=vm+o . 5* ( VXI i+vxi ) *td ( i ) 

98 CONTINUE 
AL=(XM-XT(!CI) ) /(QM*VXXI) 

TYPE 95,i\L 

IF(AL.LT.-O.l) GO TO iOO 
IF(AL.GT.O.l) GO TO 99 
GO TO 101 

99 AL=0.1 

GO TO 101 

100 AL=-0.1 

101 TYPE 105,AL,XL,YL 

105 FOPvMATC' AL=',D12.5,' XL=',D12.5, ' YL=',D12.5) 

PAUSE 

DO 120 1=2,82 
DO 110 J=2,32 

PHI(I,J)=?III(I, J)+AL*((YY(J)+YL)/((XX(I)-XL)**2 

1 +(YY(J)+YL)**2)-(YL-YY(J))/((XX(I)-:a)**2 

2 +(YL-YY(J))**2)) 

110 CONTINUE 

120 CONTINUE 
KT=0 

DT ( 1 ) =DY/ ( 1 00*VY0 (KK) ) 

121 KT=KT+1 

CALL DVOGEL(KK) 

IF(Y(KT+1) .CT.YO(KK)) GO TO 121 
SF= ( X ( KT+ 1) -XO ( KK) ) / ( XT ( KK) -XO ( KX) ) 

TYPE 122, SF 
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122 FORMAT (' SF=',D12.6) 

DO 125 1=2,82 
DO 124 J=2,32 
PHI(I,J)=SF*PHI(I, J) 

124 CONTIMUE 

125 CONTINUE 
GO TO 60 

130 DO 134 1=10,70 

TYPE 132, PHI(I,3),PHI(I,5),PHI(I,7),PHI(I,9),PHI(I,11), 
1 PHI(I, 13) ,PHI(I,15) ,PHI(I,17),PHI(I, 19) ,PHI(I,21) 

132 FORMAT(' ',10D12.5) 

134 CONTINUE 

CALL VPLOT 
GO TO 160 

140 OPEN (UNIT=1,FILE='TRAJ. DAT') 

READ (1,4) KK 
DO 150 K=1,KK 
READ(1,2) X0(K),VX0(K) 

READ(1,2) Y0(K),VY0(K) 

READ(1,2) XT(K),YT(K) 
type 155,x0(k) ,xt(k) 

155 format(' x0=',dl2.6 ,' xt=',dl2.6) 

150 CONTINUE 

GO TO 9 
160 STOP 

END 

SUBROUTINE NEWTON (X,Y, EX, EY.PSI, IK) 

IMPLICIT REAL *8(A-H,0-Z) 

COMMON/COORD/PHI(82,32) ,XX(82) ,YY(32) ,DX,DY 
DIMENSION FM(4,4,4) ,FMN(4,4,4) ,XM(4) ,XM1(4) ,YN(4) ,YN1(4) 
5 FORJIAT(2D12.6) 

IF(X.GT.XX(82)) GO TO 220 
IF(X.LT.XXd)) GO TO 220 
IF(Y.GT.YY(32)) GO TO 220 
IF(Y.LT.YYd)) GO TO 225 
HX=DX/80.0 
HY=DY/20.0 
1=0 

10 1 = 1+1 

IF(X.GT.XXd)) GO TO 10 
IF(I.LE.3) GO TO 12 
IB=I-2 
GO TO 14 
12 IB=2 

14 IF(I.EQ.82) GO TO 16 

IE=I+1 
GO TO 18 
16 IE=I 

18 J=0 

20 J=J+1 

IF(Y.GT.YY(J)) GO TO 20 
IF(J.LE.3) GO TO 25 
JB=J-2 
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GO TO 30 
25 JB=2 

30 IF(J.LE.32) GO TO 35 

JE=J+1 
GO TO 40 
35 JE=32 

40 MM=IE-IB+1 

XM(1)=1.0 
XM1(1)=0.0 
DO 70 M=2,MM 
IM=IB+M-2 

XM(M)=XM(M-1)*(X-XX(IM)) 

XM1(M)=0.0 
DO 60 1=2, M 
Xl=1.0 
DO 50 K=2,M 
IF(K.EQ.I) GO TO 50 
Xl=Xl*(X-XX(IB+K-2) ) 

50 CONTINUE 

XML(M) =XML(M)+X1 
60 CONTINUE 

70 CONTINUE 

NN=JE-JB+1 
YN(1)=1.0 
YN1(1)=0.0 
DO 100 N=2,NH 
JN=JR+N-2 

YN(N)=YN(N-1)*(Y-YY(JN)) 

YN1(N)=0.0 

DO 90 J=2,M 

Yl=1.0 

DO 80 L=2,N 

IF(L.EQ.J) GO TO 80 

Yl=Yl*(Y-YY(JB+L-2)) 

80 CONTINUE 

YN1(N)=YN1(N)+Y1 
90 CONTINUE 

100 CONTINUE 

DO 110 H=1,MM 
DO 105 N=1,NN 

FM( 1 ,M ,N ) =PH1 ( IB+M- 1 , JB+N- 1 ) 

105 CONTINUE 

no CONTINUE 

DO 140 L=2,NN 
DO 130 N=L,NN 
DO 120 M=1,MM 

FM(L,M,N) = (FM(L-1,M,N-1)-FM(L-1,M,N))/(HY*(1-D) 
1 20 CONTINUE 

130 CONTINUE 

140 CONTINUE 

DO 160 H=1,MM 
DO 150 N=1,NN 
F!-IN(1,M,N)=FM(N,M,N) 
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150 CONTINUE 

160 CONTINUE 

DO 190 K=2,MM 
DO 180 

DO 170 N=1,NN 

FMN(K,M,N)=(FMN(K-1,M-1,N)-F?IN(K-1,M,N))/(HX*(1-K)) 

170 CONTINUE 

180 CONTINUE 

190 CONTINUE 

EX=0 . 0 
EY=0.0 
PS 1=0.0 
DO 210 N=1,NN 

DO 200 M=l,l'tM 

EX=EX-XH1 (M) *YN(N)*FMN(M,M,N) 

EY=EY-YN1(N)*XM(M) *FMN(M,M,N) 

PSI=PSI+XM(M) *YN(N)*FMN(M,M,M) 

200 CONTINUE 

210 CONTINUE 

GO TO 230 
220 IK=-1 

GO TO 230 
225 IK=2 

230 RETURN 

END 

SUBROUTINE DVOGEL(K) 

IMPLICIT REAL *8(A-H,0-Z) 

COMMON/COORD/PHI (82, 32) ,XX(82) ,YY(32) ,DX,DY 
C0MM0N/DATA/X0(50) ,Y0(50) ,VX0(50) ,VY0(50) ,XT(50) ,YT(50) 
COMMON/SAVE/X( 1 10) ,Y( 1 10) ,DT( 250) ,KT 
DIM=DSQRT(DX**2+DY**2) 

QM=1.76D11 

IK=0 

IF(KT.GT.l) GO TO 200 
VX=VX0(K) 

VY=VY0(K) 

X(1)=X0(K) 

Y(1)=Y0(K) 

100 CALL NEWTON(X(KT) ,Y(KT) ,EX,EY,PSI,IK) 

101 FORl-IATC' EX=' ,D12.4,'EY=',D12.4) 

IF(IK.NE.O) GO TO 400 
AX=QM*EX 

AY=QM*EY 

XH=X(KT)-(0.5*VX-0.125*AX*DT(KT))*DT(KT) 

YH=Y(KT)-(0.5*VY-0.125*AY*DT(KT))*DT(KT) 

CALL NEWTON(XH,YH,EX,EY,PSI,IK) 

IF(IK.NE.O) GO TO 400 

AXH=QM*EX 

AYH=QM*EY 

200 XH=X(KT)+(0.5*VX+(4.0 * AX- AXH )*DT(KT)/24.0)*DT(KT) 

YH=Y(KT)+(0.5*VY+(4.0*AY-AYH)*DT(KT)/24.0)*DT(KT) 

CALL NEUTON ( XH , YH , EX , EY , PS I , IK ) 

IF(IK.NE.O) GO TO 400 
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AXH=QM*EX 

AYH=QM*EY 

X(KT+1)=X(KT)+(VX+(AX+2.0*AXH)*DT(KT)/6.0)*DT(KT) 

Y(KT+1)=Y(KT)+(VY+(AY+2.0*AYH)*DT(KT)/6.0)*DT(KT) 

260 FORMAT (21) 

AAX=AX 

AAY=AY 

CALL NEWT0N(X(KT+1) ,Y(KT+1 ) ,EX, EY, PSI , IK) 

IF(IK.NE.O) GO TO 400 

AX=QM*EX 

AY=QM*EY 

A1SQ=AAX**2+AAY**2 
IF(AlSQ.LT.l.D-6) GO TO 300 

DELD=DSQRT((X(KT+1)-X(KT))**2+(Y(KT+1)-Y(KT))**2) 

275 FORILVTC' DELD=' ,012.6/ DIM=',D12.6) 

IF(DELD.GT. DIM/40) GO TO 320 
IF(DELD.LT.DIM/120) GO TO 330 
280 VX=VX+DT(KT)*(AX+4.0*AXH+AAX)/6.0 

VY=VY+DT ( KT ) * ( AY+4 . 0*AYH+AAY) / 6 . 0 
DT(KT+1)=DT(KT) 

GO TO 440 

300 V1SQ=VX**2+VY**2 

DVSQ=AlSq*DT(KT)**2 
IF(V1SQ.GT.DVSQ*1.D4) GO TO 280 
IK=l 

GO TO 410 

320 DT(KT)=DT(KT)/1.1 

GO TO 100 

330 DT(KT)=DT(KT)*1.1 

GO TO 100 

400 IF(IK.EQ.-I) TYPE 420 

GO TO 440 
410 TYPE 430 

420 FORMAT (' ELECTRON HAS ESCAPED POTENTIAL REGION 

1 INCLUDED IN X-Y COORDINATE GRID') 

430 FORMAT (' PARTICLE VELOCITY AND ACCELERATION BOTH 

1 EQUAL ZERO') 

440 RETURN 

END 

SUBROUTINE VPLOT 
IMPLICIT REAL *8 (A-H,0-Z) 

CO>EION/COORD/PHI(82,32) ,>LX(82) ,YY(32) ,DX,DY 
DIMENSION X(250),Y(250) 

CALL INITT 
CALL SELINI 

CALL DWI NDO ( SNGL ( :iX ( 1 ) ) , SNGL ( XX ( 8 2 ) ) , SNGL ( Y Y ( 1 ) ) , SNGL ( Y Y ( 2 2 ) ) ) 
II=INT(SNGL(780*DY/DX) ) 

CALL TWniDO(0, 780,0,11) 

CALL MOVEA(SNGL(XX(2)) ,SNCL(YY(2))) 

CiVLL DRAUA( SNGL ( ( 8 2 ) ) , SNGL ( YY ( 2 ) ) ) 

PHIMIN=1.0D38 
PHrL\X=-1.0D38 
DO 10 1=2,82 
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DO 5 J=2,22 

IF(PHI(I, J) .GT.PHIMIN) GO TO 2 
PHIMIN=PHI(1,J) 

2 IF(PHI(I,J) .LT.PHIMAX) GO TO 5 

PHL'-IAX=PHI(I, J) 

5 COIITINUE 

10 CONTIIJUE 

K=0 

12 K=K+1 

V=PHIM1N+K*0 . 2 
L=0 

DO 20 J=2,22 
DO 15 1=3,82 

IF((PHI(I-1,J)-V)*(PHI(I,J)-V).GE.0.0) GO TO 15 
L=L+1 

X(L)=XX(I-1)+(XX(I)-XX(I-1))*(V-PHI(I-1,J))/(PHI(I,J)- 
1 PHI(I-l.J)) 

Y(L)=YY(J) 

15 CONTINUE 

20 CONTINUE 

DO 30 1=2,82 
DO 25 J=3,22 

IF((PHI(I,J-1)~V)*(PHI(I,J)-V).GE.O.O) GO TO 25 
L=L+1 

X(L)=XX(I) 

Y(L)=YY(J-1)+(YY(J)-YY(J-1))*(V-PHI(I,J-1))/(PHI(I,J)- 
1 PHI(I,J-D) 

25 CONTINUE 

30 CONTINUE 

LMAX=L 
XMIN=XX(82) 

DO 35 L=1,LMAX 

IF(X(L) .GT.XMIM) CO TO 35 

XMIN=X(L) 

LMIN=L 

35 CONTINUE 

XS=X(1) 

YS=Y(1) 

X(1)=X(LMIN) 

Y(1)=Y(LMIN) 

X(LMIN)=XS 

Y(U1IN)=YS 

CALL MOVEA( SNGL (X ( 1 ) ) , SNGL ( Y ( 1 ) ) ) 

LL=1 

40 LL=LL+1 

DSQ=1.0D38 
DO 45 L=LL,LMAX 

DSQL=(X(L)-X(LL-1))**2+(Y(L)-Y(LL-1))**2 

IF(DSQL.GT.DSQ) CO TO 45 

DSQ=DSQL 

LDMIN=L 

45 CONTINUE 

XS=X(LL) 



YS=Y(LL) 

X(LL)=X(LDMIN) 

Y(LL)=Y(LDMIN) 

X(LDMIN)=XS 

Y(LDMIN)=YS 

CALL DRAUA( SNGL ( X(LL) ) , SNGL ( Y (LL) ) ) 
IF(LL.EQ.LHAX) GO TO 50 
GO TO 40 

IF(V+0.2.LE.PHIMAX) GO TO 12 
CALL FINITT(800,500) 

RETURN 
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cROUTINE RTEHT.F4 

IMPLICIT REAL *8(A-H,0-Z) 

CO^GION/COORD/PHI(82,32) ,XX(82) ,YY(32) ,DX,DY 
COMMON/DATA/XO(50) ,Y0(50) ,VX0(50) ,VY0(50) ,XT(50) ,YT(50) 
COMMON/SAVE/X(110) ,Y( 110) ,DT(250) ,KT 

DIMENSION XS(llO) ,YS ( 1 10) , SX( 1 10) ,SY(110) ,TD(110) ,DDX(50) ,XM(50) 
1 ,YN(50) ,XN(50) 

EXTERNAL TEKHAN , PLTHAN , SCPHAN 
CALL INITTG 
TYPE 1 

1 FORMAT(' ENTER X AND Y DIMENSIONS OF POTENTIAL REGION', $) 

ACCEPT 2,DX,DY 

2 FORMAT (2D) 

TYPE 19 

19 FORMAT (' ENTER 1 TO INPUT DATA FROM FILE TRAJ.DAT, ENTER 0 
1 TO INPUT DATA MANUALLY' , $) 

ACCEPT 4, NT 
IF(NT.EQ.l) GO TO 140 
TYPE 3 

3 FORMJ^TC' ENTER NUl-IBER OF ELECTRON TRAJECTORIES' , $) 

ACCEPT 4,KK 

4 FORI'IAT(I) 

DO 8 K=1,KK 
TYPE 5 

5 FORILATC' ENTER INITIAL X COORDINATE AI'ID VELOCITY', $) 

ACCEPT 2,X0(K) ,VX0(K) 

TYPE 6 

6 FORMAT(' ENTER INITIAL Y COORDINATE AND VELOCITY', $) 

ACCEPT 2,Y0(K) ,VY0(K) 

TYPE 7 

7 FORI-IATC' ENTER FINAL X AIE) Y COORDINATES' , $) 

ACCEPT 2,XT(K) ,YT(K) 

type 67,k,x0(k) ,xt (k) 

67 forraat(i2,2dl2.6) 

8 CONTINUE 

9 EQ=1.60210D-19 
EM=9.1091D-31 
QM=EQ/EM 

A0=2*VX0 (KK) *VY0 (KK) / (QM* (XT (KK)-XO (KK) ) ) 

TYPE 15, AO 

15 FORI-IAT(' A0=',D) 

XX(l)=-DX/80 
XX(2)=0.0 
DO 20 1=3,82 
XX(I)=XX(I-l)+DX/80 

20 CONTINUE 
YY(l)=-DY/20 
YY(2)=0.0 

DO 30 J=3,32 
YY(J)=YY(J-l)+DY/20 
30 CONTINUE 

32 DO 50 1=1,32 

DO 40 J=l,32 
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PHI(I,J)=AO*YY(J) 

40 CONTINUE 

50 CONTINUE 

PAUSE 
IT=0 

60 DD=0.0 

D=0.0 

TYPE 65, IT 

65 F0RMAT(' ITERATION=' , 14) 

IT=IT+1 

DO 90 K=1,KK-1 

YN(K)=0.0 

KT=0 

TYPE 82,X0(K) ,Y0(K) 

DT ( 1 ) =DY/ ( 100*VY0 (K) ) 

70 KT=KT+1 

CALL DVOGEL(K) 

IF(Y(KT+1).LT.YN(K)) GO TO 75 

Y1KK)=Y(KT+1) 

aUK)=X(KT+l) 

75 IF(Y(KT+1) .GT.YO(K)) GO TO 70 

80 TYPE 82,X(KT+1) ,Y(KT+1) 

DDX(X)=(XT(K)-X(KT+1))/DABS(XT(K)-X0(K)) 

82 FORIUVT(' X=',D12.4,' Y=',D12.4) 

DD=Dn+DABS(XT(K)-X(KT+l)) 
IF(D.GT.DABS(DDX(K))) GO TO 90 
D=DABS(DDX(K)) 

XI'1(K)=X(KT+1) 

KJI=K 

DO 85 KI=1,KT+1 
XS(KI)=X(KI) 

YS(KI)=Y(KI) 

SX(KI)=FX(KI) 

SY(KI)=FY(KI) 

TD(KI)=DT(KI) 

85 CONTINUE 

90 CONTINUE 

TYPE 95,D 
TYPE 95, DD 
95 FORMAT (D) 

IF(DD.LT.0.01*DX) GO TO 130 
k=kni 

180 k=k+l 

if(k.eq.kk) go to 220 
if (ddx(k)*ddx(km) .gt .0.0) go to 180 
190 k=k+l 

if(k.eq.kk) go to 200 
if (ddx(k)*ddx(km) .gt .0.0) go to 200 
if (dabs (ddx(k) ) .It .dabs (ddx(k-l) ) ) go to 200 
go to 190 
200 lp=k-l-km 

go to 230 
lp=kk 


220 
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230 k=km 

240 k=k-l 

if(k.eq.O) go to 270 
if (ddx(k)*ddx(kia) .gt .0.0) go to 240 
250 k=k-l 

if(k.eq.O) go to 260 
if (ddx(k)*ddx(km) .gt .0.0) go to 260 
if (dabs (ddx(k) ) .It .dabs (ddx(k+l) ) ) go to 260 
go to 250 
260 ln=k+l-km 

go to 280 
270 ln=-kk 

280 if(lp.ne.kk) go to 340 

if (In.ne.-kk) go to 340 
DDN=0.0 
DO 285 K=l,KIi 
DDN=DDN+DDX(K) 

285 CONTINUE 

DDP=0.0 

DO 290 K=KM,KK-1 
DDP=DDP+DDX(K) 

290 CONTINUE 

IF(DABS(DDN) .GT.DABS(DDP)) GO TO 300 
XL=XN(1)-0.5*(XN(2)-XN(D) 

GO TO 365 

300 XL=XN(KK-l)+0.5*(XN(XK-l)-XN(KK-2)) 

GO TO 365 

340 IF(LP.GT.-LN) GO TO 350 

L=LP 

GO TO 360 
350 L=LM 

360 XL=XN(KM)+0.5*(XN(IGI+L)-XN(KM) ) 

365 YL=(1+0.1)*DY 

DDE=1.0D30 

370 YL=YL+0.05*DY 

Dl=((YN(KM)-YL)**3-3*(YN(KM)-YL)*(XN(iai)-XL)**2)/ 
1 ( ( XN ( KZi) -:a ) ** 2+ ( YN ( KM) -YL) **2 ) ** 3 

D2= ( ( Yl] ( KM) +YL ) ** 3- 3* ( YTn ( KM) +YL) * (XN ( ICi ) -XL ) ** 2 ) / 
1 ((XN(KIi)-XL)**2+(YN(KM)+YL)**2)**3 

DDA=DDE 

DDE=DABS(D1-D2) 

IF(DDE.LT.DDA) GO TO 370 

V:a 1=0.0 

DVXI 1=0.0 

VXXI=0.0 

DO 98 1=1, KT 

DVXI=DVXII 

DVXII=2*(XS(I)-XL)*((YL-YS(I))/((XS(I)-XL)**2 

1 +(YL-YS(I))**2)**2-(YS(I)+YL)/((XS(I)-XL)**2 

2 +(YS(I)+YL)**2)**2) 

VXI=VXII 

VXI I=VXI+0 . 5* (DVXI I+DVXI ) *TD ( I ) 

VXXI=VXXI+0 . 5* (VXI I+VXI ) *TD ( I ) 



98 CONTINUE 

AL= ( XII ( KM) -XT ( KM) ) / ( QM*VXXI ) 

TYPE 95, AL 

if (al . It .-0. 1) go to 100 
if (al .gt .0 . 1) go to 99 
GO TO 101 

99 AL=0.1 

GO TO 101 

100 AL=-0.1 

101 TYPE 105,AL,XL,YL 

105 F0RMAT(' AL=',D12.5,' XL=',D12.5, ' YL=',D12.5) 

DO 120 1=2,82 
DO 110 J=2,32 

Pllld, J)=PHI(I, J)+AL*((YY(J)+YL)/((XX(I)-XL)**2 

1 +(YY(J)+YL)**2)-(YL-YY(J))/((XX(I)-XL)**2 

2 +(YL-YY(J))**2)) 

110 CONTINUE 

120 CONTINUE 
KT=0 

DT(1)=DY/(100*VY0(KK)) 

121 KT=KT+1 

CALL DVOGEL(KK) 

IP(Y(KT+1) .GT.Y0(KK,)) GO TO 121 
SF=(X(KT+1)-XO(KK))/(XT(KK)-XO(KX) ) 

TYPE 122, SF 

122 FORIIATC' SF=',D12.6) 

DO 125 1=2,82 

DO 124 J=2,32 
PHI(I,J)=SF*PHI(I,J) 

124 CONTINUE 

125 CONTINUE 
GO TO 60 

130 DO 134 1=10,70 

TYPE 132, PHI(I,3),PHI(I,5),PHI(I,7),PHI(I,9),PHI(I,11), 
1 PHI(I,13) ,PHI(I,15),PHI(I,17) , PHI (I, 19) ,PHI(I,21) 

132 FOR>L\T(' ',10012. 5) 

134 CONTINUE 

CALL VPLOT 
GO TO 160 

140 OPEN (UN IT= 1 , FILE= ' TRAJ . DAT ' ) 

READ (1,4) KK 
DO 150 K=1,KK 
READ(1,2) X0(K),VX0(K) 

READ(1,2) Y0(K),VY0(K) 

READ (1,2) XT(K),YT(K) 
type 155 ,x0 (k) ,xt (k) 

155 forinat(' x0=',dl2.6 ,' xt=',dl2.6) 

150 CONTINUE 

GO TO 9 
160 STOP 

END 

SUBROUTINE NEUTON ( X , Y , EX , EY , PS I , IK) 

IMPLICIT REAL *8(A-H,0-Z) 



COIC'ION/ COORD / PHI(82,32),XX(82),YY(32),DX,DY 
DIMENSION FM(4,4,4) ,FMN(4,4,4) ,XM(4) ,XJ-I1(4) ,YN(4) ,YN1(4) 
5 FORI»!AT(2D12.6) 

IF(X.GT.XX(82)) GO TO 220 
IF(X.LT.XXd)) GO TO 220 
IF(Y.GT.YY02) ) GO TO 220 
IF(Y.LT.YYd)) GO TO 225 
HX=DX/80.0 
HY=DY/20.0 
1=0 

10 1 = 1+1 

IF(X.GT.XXd)) GO TO 10 
IF(I.LE.3) GO TO 12 
IB=I-2 
GO TO 14 
12 IB=2 

14 IF(I.EQ.82) GO TO 16 

IE=I+1 
GO TO 18 
16 IE=I 

18 J=0 

20 J=J+1 

IF(Y.GT.YY(J)) GO TO 20 
IF(J.LE.3) GO TO 25 
JB=J-2 
GO TO 30 
25 JB=2 

30 IF(J.EQ.32) GO TO 35 

JE=J+1 
GO TO 40 
35 JE=J 

40 MM=IE-IB+1 

XM(1)=1.0 
XM1(1)=0.0 
DO 70 M=2,MM 
IM=IB+M-2 

XI‘l(M)=XM(M-l)*(X-XXdN) ) 

XM1(M)=0.0 
DO 60 1=2 ,M 
Xl=1.0 
DO 50 K=2,M 
IF(K.EQ.I) GO TO 50 
X 1=X 1 * ( X-XX ( IB+K- 2 ) ) 

50 CONTINUE 

XM1(M)=XM1(M)+X1 
60 CONTINUE 

70 CONTINUE 

NN=JE-JB+1 
YN(1)=1.0 
YN1(1)=0.0 
DO 100 N=2,NN 
JN=JB+N-2 

YM (N ) =YN ( N- 1 ) * ( Y- YY ( JN ) ) 
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YN1(N)=0.0 

DO 90 J=2,N 

Yl=1.0 

DO 80 L=2,N 

IF(L.EQ.J) GO TO 80 

Yl=Yl*(Y-YY(JB+L-2)) 

80 CONTINUE 

YN1(N)=YN1(N)+Y1 
90 CONTINUE 

100 CONTINUE 

DO 110 M=1,MM 
DO 105 N=1,NN 

FM ( 1 , M , N ) =PHI ( IB+M- 1 , JB+M- 1 ) 

105 CONTINUE 

110 CONTINUE 

DO 140 L=2,NN 
DO 130 N=L,NN 
DO 120 M=1,MM 

FMa,M,N) = (FM(L-l ,M,N-1 )-FM(L-l ,M,M) ) / (HY*( 1-L) ) 

120 CONTIinJE 

130 CONTINUE 

140 CONTINUE 

DO 160 M=1,NM 
DO 150 N=1,NM 
FNN(l,M,N)=FM(N,n,N) 

150 CONTINUE 

160 CONTINUE 

DO 190 K=2,MM 

DO 180 M=K,MM 

DO 170 N=1,NN 

FMN ( K , M , N ) = ( Fl-DI ( K- 1 , M- 1 , N ) - F? m ( K- 1 , M , N ) ) / ( HX * ( 1 -K ) ) 
170 CONTINUE 

180 CONTINUE 

190 CONTINUE 

EX=0 . 0 
EY=0.0 
PSI=0.0 
DO 210 N=1,NN 

DO 200 M=1,MM 

EX=EX-XM1 (M) *Y1UN) *FNN(M,M,N) 

EY=EY-YN1 (N) *XH(M) *FMN(M,M,N) 

PS I=PS I+X!1(M) *YN ( N ) *FMJI (M ,M , N ) 

200 CONTINUE 

210 CONTINUE 

GO TO 230 
220 IK=-1 

GO TO 230 
225 IK=2 

230 RETURN 

END 

SUBROUTINE DVOGEL(K) 

IMPLICIT REAL *8(A-H,0-Z) 

CO?!MON/COORD/PHI(82,32) ,XX(82) ,YY(32) ,DX,DY 
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COJ^MON/DATA/XO(50) ,Y0 (50) ,VX0(50) ,VY0(50) ,XT(50) ,YT(50) 
COMHON/SAVE/X(110) ,Y(110) ,DT(250) ,KT 
DIM=DSQRT(DX**2+DY**2) 

QM=1.76D11 

IK=0 

IF(KT.GT.l) GO TO 200 
VX=VX0(K) 

VY=VY0(K) 

X(1)=X0(K) 

Y(1)=Y0(K) 

100 CALL NEOTON(X(KT) ,Y(KT),EX,EY,PSI,I!0 

101 FOR!tAT(' EX=',D12.4/EY=',D12.4) 

IF(IK.NE.O) GO TO 400 
AX=QM*EX 

AY=QM*EY 

XH=X(KT)-(0.5*VX-0.125*AX*DT(KT))*DT(KT) 

YH=Y(XT)-(0.5*VY-0.125*AY*DT(KT))*DT(KT) 

CALL NEWT0N(XH,YH,EX,EY,PSI,IK) 

IF(IK.NE.O) GO TO 400 

AXH=QM*EX 

AYH=QM*EY 

200 Iffi=X ( KT ) + ( 0 . 5 * VX+ ( 4 . 0* AX-AXH) *DT ( KT ) / 2 4 . 0 ) *DT ( KT ) 

YH=Y(KT)+(0.5*VY+(4.0*AY-AYH)*DT(KT)/24.0)*DT(KT) 

CALL NEWTON(XH,YH,EX,EY,PSI,IK) 

IF(IK.NE.O) GO TO 400 

AXH=QM*EX 

AYH=QM*EY 

X(KT+l)=X(KT) + (VX+(AX+2.0*AXIl)*DT(KT)/6.0)*DT(KT) 

Y ( KT+ 1 ) = Y ( KT ) + ( VY+ ( AY+2 . 0*AYH ) *DT ( KT ) / 6 . 0 ) *DT ( KT ) 

260 F0RHAT(2I) 

AAX=AX 

AAY=AY 

CALL NEVJTOM(X(KT+l) ,YaCT+l) ,EX,EY,PSI,IK) 

IF(IK.ME.O) GO TO 400 

AX=QH*EX 

AY=QM*EY 

A1SQ=AAX**2+AAY**2 
IF(AlSQ.LT.l.D-6) GO TO 300 

DELD=D SQRT ( ( X ( KT+1 ) -X ( KT) ) ** 2+ ( Y ( KT+1 ) -Y ( KT ) ) ** 2 ) 

275 FORIIATC' DELD=' ,012.6/ DIM=',D12.6) 

IF(DELD.GT. DIM/40) GO TO 320 
IF(DELD.LT. DIM/120) CO TO 330 
280 VX=VX+DT ( KT ) * ( AX+4 . 0* AXH+Ai\X) / 6 . 0 

VY= VY+DT ( KT ) * ( AY+4 . 0* AYH+AAY ) / 6 . 0 
DT(KT+1)=DT(KT) 

GO TO 440 

300 V1SQ=VX**2+VY**2 

DVSQ=A1SQ*DT(KT)**2 
IF(VlSQ.GT.DVSq*l.D4) CO TO 280 
IK=1 

GO TO 410 
DT(KT)=DT(KT) /I. 1 
GO TO 100 


320 
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330 DT(KT)=DT(KT)*1 .1 

GO TO 100 

400 IF(IK.EQ.-l) TYPE 420 

GO TO 440 
410 TYPE 430 

420 FORMAT(' ELECTRON HAS ESCAPED POTENTIAL REGION 

1 INCLUDED IN X-Y COORDINATE GRID') 

430 FORMAT (' PARTICLE VELOCITY AND ACCELERATION BOTH 

1 EQUAL ZERO') 

440 RETURN 

END 

SUBROUTINE VPLOT 
IMPLICIT REAL *8 (A-H.O-Z) 
COMMON/COORD/PHI(82,32),XX(82),YY(32),DX,DY 
DIMENSION X(250),Y(250) 

CALL INITT 
CALL SELINI 

CALL DWIi'fDO(SNGL(XX(l)),SNGL(XX(82)) ,SNGL(YY(D) ,SNGL(YY(22))) 
II=INT ( SNGL ( 780*DY/DX) ) 

CALL TWINDO(0, 780, 0,11) 

CALL MO VEA( SNGL ( XX( 2 ) ) , SNGL ( YY ( 2 ) ) ) 

CALL DRAUA ( SNGL ( XX ( 8 2 ) ) , SNGL ( YY ( 2 ) ) ) 

PHIMIN=1.0D38 
PHIMAX=-1.0D38 
DO 10 1=2,82 
DO 5 J=2,22 

IF(PIII(I,J) .GT.PHIMIN) GO TO 2 
PHIMIN=PHI(I,J) 

2 IF(PHI(I,J) .LT,PHIMi\X) CO TO 5 
PHIMAX=PHI(I,J) 

5 CONTIJTUE 
10 CONTINUE 
K=0 

12 K=K+1 

V=PHIMIN+K*0 . 2 
L=0 

DO 20 J=2,22 
DO 15 1=3,32 

IF((PHI(I-1,J)-V)*(PHI(I,J)-V).GE.0.0) GO TO 15 
L=L+1 

X(L)=XX(I-1)+(XX(I)-XX(I-1))*(V-PHI(I-1,J))/(PHI(I,J)- 
1 ?HI(I-1,J)) 

Y(L)=YY(J) 

15 CONTINUE 
20 CONTINUE 

DO 30 1=2,82 
DO 25 J=3,22 

IF((PHI(I,J-1)-V)*(PHI(I,J)-V) .GE.0.0) GO TO 25 

L=L+1 

X(L)=XX(I) 

Y(L)=YY(J-1) + (YY(J)-YY(J-1))*(V-PHI(I,J-1) )/(PlII(I, J)- 
1 PHI(I,J-D) 

CONTINUE 


25 
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30 COIITINUE 
LMAX=L 
X?'IIN=XX(82) 

DO 35 L=1,LMAX 

IF(X(L) .GT.XHIN) GO TO 35 

XMIN=X(L) 

UTIN=L 

35 CONTINUE 

XS=X(1) 

YS=Y(1) 

X(1)=X(LMIN) 

Y(1)=Y(LMII1) 

X(LMIN)=XS 

Y(LMIN)=YS 

CALL MOVEA( SNGL ( X (1 ) ) , SNGL (Y ( 1) ) ) 

LL=1 

40 LL=LL+1 

DSQ=1.0D38 
DO 45 L=LL,LMAX 

DSQL=(X(L)-X(LL-1))**2+(Y(L)-Y(LL-1))**2 

IF(DSQL.GT.DSQ) GO TO 45 

DSQ=DSQL 

LDHIN=L 

45 CONTINUE ' 

XS=X(LL) 

YS=Y(LL) 

X(LL)=X(LDMIN) 

Y(LL)=Y(LDMIN) 

X(LDMIN)=XS 

Y(LDHIN)=YS 

CALL DRAUA(SNGL(X(LL) ) ,SNGL(Y(LL) ) ) 
IF(LL.EQ.LMAX) GO TO 50 
GO TO 40 

50 IF(V+0.2.LE.PHIHAX) GO TO 12 

CALL FIMITT(800,500) 

RETUFJT 

E?ro 
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APPENDIX E 

Fortran Routines for Quadrupole Synthesis of 
Solutions to The Model Problem 
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cROUTINE QP0LE.F4 

IMPLICIT REAL *8(A-H,0-Z) 

COMMON/COORD/PHI(82,32) ,XX(82) ,YY(32) ,DX,DY 
COMMON/DATA/XO(50) ,Y0(50) ,VX0(50) ,VY0(50) ,XT(50) ,YT(50) 
COMMON/SAVE/X(110) ,Y(110) ,DT(250) ,KT 
DIMENSION XS(llO) ,YS(110) ,TD( 1 10) ,DDX(50) 

1 ,YN(50) ,XN(50) 

EXTERNAL TEKHAN,PLTHAN,SCPHAM 
CALL INITTG 
TYPE 1 

1 FORI-IATC' ENTER X AND Y DIMENSIONS OF POTENTIAL REGION', $) 
ACCEPT 2,DX,DY 

2 FORMAT (2D) 

TYPE 19 

19 FORMATC' ENTER 1 TO INPUT DATA FROM FILE TRAJ.DAT, ENTER 0 
1 TO INPUT DATA MANUALLY', $) 

ACCEPT 4, NT 
IF(NT.EQ.l) GO TO 140 
TYPE 3 

3 FORMATC' ENTER NUMBER OF ELECTRON TRAJECTORIES' , $) 

ACCEPT 4,KK 

4 FORllAT(I) 

DO 8 K=1,KK 
TYPE 5 

5 FORMAT(' ENTER INITIAL X COORDINATE AND VELOCITY', $) 

ACCEPT 2,X0(K) ,VXO(K) 

TYPE 6 

6 FORMATC' ENTER INITIAL Y COORDINATE AND VELOCITY',?) 

ACCEPT 2,Y0(K) ,VYOCK) 

TYPE 7 

7 FORMATC' ENTER FINAL X AND Y COORDINATES',?) 

ACCEPT 2,XT(K) ,YT(K) 

type 67,k,xO(k) ,xt(k) 

67 format (12, 2dl 2. 6) 

8 CONTINUE 

9 EQ=1.60210D-19 
EM=9.1091D-31 
QM=EQ/EM 

A0=2AVX0(KK)*VY0(KK) /(QM*(XT(KK)-X0(iaO ) ) 

TYPE 15, AO 

15 FORI-IATC' AO=',D) 

XX(l)=-DX/80 
XX(2)=0.0 
DO 20 1=3,32 
XX(I)=XX(I-l)+DX/80 

20 CONTINUE 
YY(l)=-DY/20 
YY(2)=0.0 

DO 30 J=3,32 
YY(J)=YY(J-l)+DY/20 
30 CONTINUE 

32 DO 50 1=1,32 

DO 40 J=l,32 
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PHI(I,J)=AO*YY(J) 

40 CONTITTOE 

50 CONTINUE 

PAUSE 
IT=0 

60 DD=0.0 

D=0.0 

TYPE 65, IT 

65 F0RMAT(' ITERATION=' ,14) 

IT=IT+1 

DO 90 K=1,KE-1 

YTI(K)=0.0 

KT=0 

TYPE 82,X0(K) ,Y0(K) 

DT ( 1) =DY/ ( 1 00* VYO ( K) ) 

70 KT=KT+1 

CALL DVOGEL(K) 

IF(Y(KT+1).LT.YT'I(K)) GO TO 75 
Y!J(K)=Y(KT+1) 

XN(K)=X(KT+1) 

75 IF(Y(KT+1) .GT.YO(K)) GO TO 70 

80 TYPE 82,X(KT+1) ,Y(KT+1) 

DDX(K)=(XT(K)-X(KT+1))/(XT(K)-X0(K)) 

82 F0R1-1AT(' X=',D12.4,' Y=',D12.4) 

DD=DD+DABS ( XT ( K) -X ( KT+1 ) ) 
IF(D.GT.DABS(DDX(K))) GO TO 90 
D=DABS(DDX(K)) 

XM=X(KT+1) 

KM=K 

DO 85 KI=1,ET+1 
XS(KI)=X(KI) 

YS(KI)=Y(KI) 

TD(KI)=DT(KI) 

85 CONTINUE 
90 CONTINUE 
TYPE 95, D 
TYPE 95, DD 
95 FORMAT(D) 

IF(DD.LT.0.01*DX) GO TO 126 
XL=X1I(KM) 

YL=i.4*DY 
VXI 1=0.0 
DVXI 1=0.0 
VXXI=0.0 
DO 98 1=1, KT 
DVXI=DVXII 

DVXII=( (YS(I)+YL)**3-3*(YS(I)+YL)*(XS(I)-:a) 

1 **2)/((XS(I)-XL)**2+(YS(I)+YL)**2)**3 

2 +((YS(I)-YL)**3-3*(YS(I)-YL)*(XS(I)-XL)**2 

3 )/((XS(I)-:a)**2+(YS(I)-YL)**2)**3 
v:\i=vxi I 

VXII=VXI+0 . 5* (DVXI I+DVXI) *TD ( I) 

V>LXI = VXXI+0 . 5 * ( VXI I+VXI ) *TD ( I ) 
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98 CONTINUE 

AL= ( XM-XT ( 101) ) / ( QM* VXXI ) 

TYPE 95, AL 

IF(AL.LT.-O.l) GO TO 100 
IF(AL.GT.O.l) GO TO 99. 

GO TO 101 

99 AL=0.1 

GO TO 101 

100 AL=-0.1 

101 TYPE 105,AL,XL,YL 

105 FORIIATC' AL=',D12.5,' XL=',D12.5, ' YL=',D12.5) 

DO 120 1=2,82 
DO 110 J=2,32 

PHI(I,J)=PHI(I,J)+AL*(XX(I)-XL)*((YY(J)+YL)/(( 

1 XX(I)-XL)**2+(YY(J)+YL)**2)**2+(YY(J)-YL)/ 

2 ((XX(I)-XL)**2+(YY(J)-YL)**2)**2) 

no CONTINUE 

120 CONTINUE 
KT=0 

DT ( 1 ) =DY/ ( 100*VY0 (KK) ) 

TYPE 82,X0(KX) ,Y0(KK) 

121 KT=KT+1 

CALL DVOGEL(KK) 

IF(Y(KT+1) .GT.YO(KK)) GO TO 121 
TYPE 82,X(KT+1) ,Y(KT+1) 

SF= ( X ( KT+1) -XO (KK) ) / (XT (XK) -XO (KK) ) 

TYPE 122,SF 

122 F0R1L\T(' SF=',D12.6) 

DO 125 1=2,82 

DO 124 J=2,32 
PHI(I,J)=SF*PHI(I,J) 

124 CONTINUE 

125 CONTINUE 
GO TO 60 

126 DO 128 1=10,70 

TYPE 129 ,PHI(I,3),PHI(I,5),PHI(I,7),PHI(I,9),PHI(I,11) 
1 ,PHI(I, 13),PHI(I,15),PHI(I,17),PHI(I,19) ,PHI(I,21) 

128 CONTINUE 

129 FORMAT (lOD 12. 4) 

130 CALL VPLOT 
CO TO 160 

140 OPEN (UMIT= 1 , FILE=' TRAJ .DAT' ) 

READ (1,4) KK 
DO 150 K=1,KK 
READ(1,2) X0(K),VX0(K) 

READ(1,2) Y0(K),VY0(K) 

READ(1,2) XT(K),YT(K) 
type 155 ,xO(k) ,xt (k) 

155 forraat(' x0=',dl2.6 ,' xt=',dl2.6) 

150 CONTINUE 

GO TO 9 
160 STOP 

END 
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SUBROUTINE NEOTON ( X , Y , EX , EY , PS I , IK) 

IMPLICIT REAL *8(A-H,0-Z) 

COMMON/COORD/PHI (82, 32) ,XX(82) ,YY(32) ,DX,DY 
DIMENSION FM(4,4,4) ,FMN(4,4,4) ,XM(4) ,XM1(4) ,YN(4) ,YN1(4) 
5 FORMAT(2D12.6) 

IF(X.GT.XX(82)) GO TO 220 
IF(X.LT.XXd)) GO TO 220 
IF(Y.GT.YY(22)) GO TO 220 
IF(Y.LT.YYd)) GO TO 225 
HX=DX/80.0 
HY=DY/20.0 
1=0 

10 1 = 1+1 

IF(X.GT.XXd)) GO TO 10 
IF(I.LE.3) GO TO 12 
IB=I-2 
GO TO 14 
12 IB=2 

14 IF(I.EQ.82) GO TO 16 

IE=I+1 
GO TO 18 
16 IE=I 

18 J=0 

20 J=J+1 

IF(Y.GT.YY(J)) GO TO 20 
IF(J.LE.3) GO TO 25 
JB=J-2 
GO TO 30 
25 JB=2 

30 IF(J.EQ.32) GO TO 35 

JE=J+1 
GO TO 40 
35 JE=J 

40 MM=IE-IB+1 

XM(1)=1.0 
XM1(1)=0.0 
DO 70 M=2,MM 
IM=IB+M-2 

XM (M) =XM(M- 1 ) * ( X-XX ( IM) ) 

XM1(M)=0.0 

DO 60 1=2 ,M 

Xl=1.0 

DO 50 K=2,M 

IF(K.EQ.I) CO TO 50 

Xl=Xl*(X-XX(IB+K-2)) 

50 CONTINUE 

XM1(M)=XM1(M)+X1 
60 CONTINUE 

70 CONTINUE 

NN=JE-JR+1 
YN(1)=1.0 
YN1(1)=0.0 
DO 100 N=2,NN 



JN=JB+N-2 

TO ( N ) =YN ( N- 1 ) * ( Y-Y Y ( Jl'i ) ) 

TOl(N)=0.n 
DO 90 J=2,N 
Yl=1.0 
DO 80 L=2,N 
IF(L.EQ.J) GO TO 80 
Yl=Yl*(Y-YY(JB+L-2) ) 

80 CONTINUE 

TO1(N)=TO1(N)+Y1 
90 CONTINUE 

100 CONTINUE 

DO 110 M=1,MM 
DO 105 N=1,NN 

EM (1 ,M , N ) =PHI ( IB+M- 1 , JB+M- 1 ) 

105 CONTINUE 

110 CONTINUE 

DO lAO L=2,NN 
DO 130 N=L,NN 
DO 120 M=l,>tM 

FM(L,M,N) = (FM(L-l,M,N-l)-FMa-l,M,M))/(HY*(l-L)) 

120 CONTINUE 

130 CONTINUE 

140 CONTINUE 

DO 160 M=1,MM 
DO 150 N=1,NN 
FtIN(l,M,M)=FM(N,M,N) 

150 CONTINUE 

160 CONTINUE 

DO 190 K=2,;-tM 

DO 180 M=K,MI 

DO 170 N=1,NN 

FMN(K,M,N)=(FMN(K-1 ,M-1 ,N)-FNN(K- 1 ,M,N) ) / (HX* ( 1-K) ) 
170 CONTINUE 

180 CONTINUE 

190 CONTINUE 

EX=0.0 
EY=0 . 0 
PSI=0.0 
DO 210 N=1,NN 

DO 200 H=l,fDI 

EX=EX-XM1 (M) *TO(N) *FMN(M,M,N) 

EY = E Y- YN 1 ( N ) * XM ( M ) * FMN ( M , H , N ) 
PSI=PSI+XM(M)*YN(N)*FMN(M,M,N) 

200 CONTINUE 

210 CONTINUE 

CO TO 230 
220 IK=-1 

GO TO 230 
225 IK=2 

230 RETURN 

END 

SUBROUTINE D VOGEL (K) 
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IMPLICIT REAL *8(A-H,0-Z) 

CO?IMON/COORD/PHI(82,32) ,XX(82) ,YY(32) ,DX,DY 
COMMON/DATA/XO(50) ,Y0(50) ,VX0(50) ,VY0(50) ,XT(50) ,YT(50) 
COMMON/SAVE/X(110) ,Y(110) ,DT(250) ,KT 
DIM=DSQRT(DX**2+DY**2) 

QM=1.76D11 

IK=0 

IF(KT.GT.l) GO TO 200 
VX=VX0(K) 

VY=VY0(K) 

X(1)=X0(K) 

Y(1)=Y0(K) 

100 CALL NEUTON(X(KT) ,Y(KT) ,EX,EY,PSI,IK) 

101 FORMATC' EX=',D12.4/EY=',D12.4) 

IF(IK.NE.O) GO TO 400 
AX=QM*EX 

AY=QM*EY 

XH=X ( KT ) - ( 0 . 5 * VX- 0 . 1 2 5 * AX*DT ( KT ) ) *D T ( KT ) 

YH=Y ( KT ) - ( 0 . 5* VY- 0 . 1 25* AY*DT ( KT ) ) *DT ( KT ) 

CALL NEUTON(XH,YH,EX,EY,PSI,IK) 

IF(IK.NE.O) GO TO 400 

AXH=QM*EX 

AYH=QM*EY 

200 XII=X(KT)+(0.5*VX+(4.0*AX-AXH)*DT(KT)/24.0)*DT(KT) 

YH=Y(KT)+(0.5*VY+(4.0*AY-AYH)*DT(KT)/24.0)*DT(KT) 

CALL NEWTON (XH,YH, EX, EY,PSI, IK) 

IF(IK.NE.O) GO TO 400 

AXH=QM*EX 

AYH=QM*EY 

X(KT+1)=X(KT)+(VX+(AX+2.0*AXH)*DT(KT)/6.0)*DT(KT) 

Y(KT+1)=Y(KT)+(VY+(AY+2.0*AYH)*DT(KT)/6.0)*DT(KT) 

260 FORMAT(2I) 

AAX=AX 

AAY=AY 

CALL MEl>rrON(X(KT+l) ,Y(KT+1),EX,EY,PSI,IK) 

IF(IK.NE.O) GO TO 400 

AX=QM*EX 

AY=QM*EY 

A1SQ=AAX**2+AAY**2 
IF(A1SQ.LT. l.D-6) GO TO 300 

DELD=DSQRT((X(KT+1)-X(KT))**2+(Y(KT+1)-Y(KT))**2) 

275 F0R1UT(' DELD=' ,D12.6,' DIM=',D12.6) 

IF(DELD.GT. DIM/40) GO TO 320 
IF(DELD.LT. DIM/120) GO TO 330 
280 VX=VX+DT(KT)*(AX+4.0*AXH+AAX)/6.0 

VY=VY+DT ( KT ) * ( AY+4 . 0*AYH+AAY) / 6 . 0 
DT(KT+1)=DT(KT) 

GO TO 440 

300 V1SQ=VX**2+VY**2 

DVSQ=A1SQ*DT(KT)**2 
IF(V1SQ.GT.DVSQ*1.D4) GO TO 280 
IK=1 

GO TO 410 
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320 DT(KT)=DT(KT)/1.1 

GO TO 100 

330 DT(KT)=DT(KT)*1.1 

GO TO 100 

400 IF(IK.EQ.-l) TYPE 420 

GO TO 440 
410 TYPE 430 

420 FORMAT(' ELECTRON HAS ESCAPED POTENTIAL REGION 

1 INCLUDED IN X-Y COORDINATE GRID') 

430 FORMAT(' PARTICLE VELOCITY AND ACCELERATION BOTH 

1 EQUAL ZERO') 

440 RETURN 

END 

SUBROUTINE VPLOT 
IMPLICIT REAL *8 (A-H.O-Z) 
COMMON/COORD/PHI(82,32),XX(82),YY(32) ,DX,DY 
DIMENSION X(250) ,Y(250) 

CALL IN ITT 
CALL SELINI 

CALL DUINDO(SNGL(XX(l)),SNGL(XX(82)) ,SNGL(YY(D) ,SNGL(YY(22))) 
II=INT ( SNGL ( 780*DY/DX) ) 

CALL TWINDO (0,780,0,11) 

CALL MOVE A( SNGL ( XX ( 2 ) ) , SNGL ( YY ( 2 ) ) ) 

CALL DRAWA(SNGL(XX(82)) ,SNGL(YY(2))) 

PHIMIN=1.0D38 
PHIMAX=-1.0D38 
DO 10 1=2,82 
DO 5 J=2,22 

1F(PHI(I,J) .GT.PHIMIN) GO TO 2 
PH1M1N=PHI(I,J) 

2 IF(PHI(I,J).LT.PHIMAX) GO TO 5 

PHIMAX=PHI(I,J) 

5 CONTINUE 

10 CONTINUE 

K=0 

12 K=K+1 

V=PHIMIN+K*0.2 

L=0 

DO 20 J=2,22 
DO 15 1=3,82 

IF((PHI(I-1,J)-V)*(PHI(I,J)-V).GE.0.0) GO TO 15 
L=L+1 

X(L)=XX(I-1)+(XX(I)-XX(I-1))*(V-PHI(1-1,J))/(PHI(1,J)- 
1 PHI(I-1,J)) 

Y(L)=YY(J) 

15 CONTINUE 

20 CONTINUE 

DO 30 1=2,82 
DO 25 J=3,22 

IF((PHI(I,J-1)-V)*(PHI(I,J)-V).GE.0.0) GO TO 25 

L=L+1 

X(L)=XX(I) 

Y(L)=YY(J-1)+(YY(J)-YY(J-1))*(V-PHI(I,J-1))/(PHI(I,J)- 
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1 PHI(I,J-D) 

25 CONTINUE 

30 CONTINUE 

LMAX=L 
XMIN=XX(82) 

DO 35 L=1,LMAX 

IF(X(L) .GT.XMIN) GO TO 35 

XMIN=X(L) 

LMIN=L 

35 CONTINUE 

XS=X(1) 

YS=Y(1) 

X(1)=X(LMIN) 

Y(1)=Y(LMIN) 

X(LMIN)=XS 

Y(LMIN)=YS 

CALL MOVEA(SNGL(X(D) ,SNGL(Y(1))) 

LL=1 

40 LL=LL+1 

DSQ=1.0D38 
DO 45 L=LL,LMAX 

DSQL=(X(L)-X(LL-1))**2+(Y(L)-Y(LL-1))**2 

IF(DSQL.GT.DSQ) GO TO 45 

DSQ=DSQL 

LDMIN=L 

45 CONTINUE 

XS=X(LL) 

YS=Y(LL) 

X(LL)=X(LDMIN) 

Y(LL)=Y(LDMIN) 

X(LDMIN)=XS 

Y(LDMIN)=YS 

CALL DRAUA ( SNGL ( X ( LL ) ) , SNGL ( Y ( LL) ) ) 
IF(LL.EQ.LMAX) GO TO 50 
GO TO 40 

50 IF(V+0.2.LE.PHIMAX) GO TO 12 

CALL FINITT(800,500) 

RETURN 

END 
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cROUTINE QPOLES.F4 

IMPLICIT REAL *8(A-H,0-Z) 

COMfTON/COORD/PHI(82,32) ,XX(82) ,YY(32) ,DX,DY 
COMMON/DATA/XO(50) ,Y0(50) ,VX0(50) ,VY0(50) ,XT(50) ,YT(50) 
COMMON/SAVE/X(110) ,Y(110) ,DT(110) ,KI(50) 
COMMON/COEFF/VI(5,6) ,AL(5) 

DIMENSION XN(6) 

EXTERJIAL TEK!IAN,PLTHAN,SCPnAN 
CALL INITTG 
TYPE 1 

1 FORJIAT(' ENTER X AND Y DIMENSIONS OF POTENTIAL REGION', $) 
ACCEPT 2,DX,DY 

2 FORMAT (2D) 

TYPE 19 

19 FORMAT (' ENTER 1 TO INPUT DATA FROM FILE TRAJ.DAT, ENTER 0 
1 TO INPUT DATA MANUALLY', $) 

ACCEPT 4, NT 
IF(NT.EQ.l) GO TO 140 
TYPE 3 

3 FORMAT(' ENTER NUMBER OF ELECTRON TRAJECTORIES' , $) 

ACCEPT 4,KK 

4 FORMAT(I) 

DO 8 K=1,KK 
TYPE 5 

5 FORMAT (' ENTER INITIAL X COORDINATE AND VELOCITY', $) 

ACCEPT 2,X0(K) ,VX0(K) 

TYPE 6 

6 FORMAK' ENTER INITIAL Y COORDINATE AND VELOCITY', $) 

ACCEPT 2,Y0(K) ,VY0(K) 

TYPE 7 

7 FORMAK' ENTER FINAL X AND Y COORDINATES' , $) 

ACCEPT 2,XT(K) ,YT(K) 

type 670,k,x0(k) ,xt (k) 

670 fornat (12 , 2dl2 . 6) 

8 CONTINUE 

9 EQ=1.60210D-19 
EM=9.1091D-31 
QM=EQ/EM 

A0=2*VX0 (KK) *VY0 (KK) / (QM* (XT (KK) -XO (KK) ) ) 

TYPE 15, AO 

15 FORMAK' A0=',D) 

XX(l)=-DX/80 
XX(2)=0.0 
DO 20 1=3,82 
XX(I)=XX(I-l)+DX/80 

20 CONTINUE 
YY(l)=-DY/20 
YY(2)=0.0 

DO 30 J=3,32 
YY(J)=YY(J-1)+DY/20 
30 CONTINUE 

32 DO 50 1=1,82 

DO 40 J=l,32 
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PHI(I,J)=AO*YY(J) 

40 CONTINUE 

50 CONTINUE 

IT=0 
ITT=0 
DD=1.0D30 

69 D=DD 
DD=0.0 
RD=0.0 

DO 90 K=1,KK 
KT=0 

YMAX=0.0 

DT(1)=DY/(100*VY0(K)) 

TYPE 80,X0(K) ,Y0(K) 

70 lCr=KT+l 

CALL DVOGEL(K,KT,KIK) 

IF(Y(KT+1) .LT.YMAX) GO TO 75 
YMAX=Y(KT+1) 

XN(K)=X(KT+1) 

75 IF(Y(KT+1) .GT.YO(K)) GO TO 70 

TYPE 80, X(KT+1),Y(KT+1) 

IF(K.EQ.KK) GO TO 90 

80 F0RMAT(' X=',D12.6,' Y=',D12.6) 

DD=DD+DABS ( XT ( K) -X(KT+1 ) ) 
RD=RD+(X(KT+1)-X0(K))/(XT(K)-X0(K)) 

85 F0RIIAT(' IT=',I3,' DD=',D12.6,' RD=',D12.6) 

YL=1.4*DY 
DO 87 L=1,KK-1 
XL=XN(L) 

VI 1=0.0 
DVI 1=0.0 
VI(K,L)=0.0 
DO 86 1=1, KT 
DVI=DVII 

DVII=((Y(I)+YL)**3-3*(Y(I)+YL)*(X(I)-XL)**2)/((X(I)-XL)**2 

1 +(Y(I)+YL)**2)**3+((Y(I)-YL)**3-3*(Y(I)-YL)*(X(I) 

2 -XL)**2)/((X(I)-XL)**2+(Y(I)-YL)**2)**3 
VVI=VII 

VI I=VVI+0 . 5* (D VI I+D VI ) *DT ( I ) 
VI(K,L)=VI(K,L)+0.5*(VII+VVI)*DT(I) 

86 CONTINUE 
VI(K,L)=VI(K,L)*QM 
TYPE 88,K,L,VI(K,L) 

87 CONTINUE 

88 FORMAT!' K=',I3,' L=',I3,' VI(K,L)=' ,D12.6) 

VI(K,KK)=X(KT+1)-XT(K) 

TYPE 89,K,VI(K,KK) 

89 FORl-lAT!' K=',I3,' DISCREPANCY=' ,D12.6) 

90 CONTINUE 

TYPE 85,IT,DD,RD 
IF(D.LT.DD) ITT=ITT+1 
IF(ITT.GT.3) GO TO 125 
IF(DD.LT. .01*DX) GO TO 126 
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ERR=D AB S ( ( X ( KT+ 1 ) -XT ( KK) ) /XT ( KK) ) 

IF(ERR.LT.0.03) GO TO 905 
SF=(X(KT+1)-X0(KK))/(XT(KK)-X0(KK)) 

TYPE 901.SF 

901 F0RMAT(' SCALE FACT0R=' ,D12.6) 

DO 903 1=2,82 

DO 902 J=2,32 
PHI(I,J)=SF*PHI(I,J) 

902 CONTINUE 

903 CONTINUE 
GO TO 69 

905 CALL GAUSS (KK-1) 

AA= .0 

DO 906 K=1,KK-1 

IF(DABS(AL(K)) .LT.AA) GO TO 906 
AA=DABS(AL(K)) 

906 CONTINUE 

DO 908 K=1,KK-1 
IF(AA.LT.O.l) GO TO 908 
AL(K)=0.1*AL(K) /AA 
908 CONTINUE 

IT=IT+1 
YL=1.4*DY 
DO 95 L=1,KK-1 
XL=XN(L) 

TYPE 91,AL(L) ,XL,YL 

91 FORMATC' AL=',D12.6,' XL=%D12.6,' YL=',D12.6) 

DO 93 1=2,82 

DO 92 J=2,32 

PHI(I, J)=PHI(I, J)+AL(L)*(XX(I)-XL)*((YY(J)+YL)/((XX(I)-XL)**2 

1 +(YY(J)+YL)**2)**2+(YY(J)-YL)/((XX(I)-XL)**2 

2 +(YY(J)-YL)**2)**2) 

92 CONTINUE 

93 CONTINUE 

95 CONTINUE 

GO TO 69 

124 FORMAT(' FOUR DEPARTURES FROM MONOTONE DECREASING 
1 SUM OF ABSOLUTE DISCREPANCIES HAVE OCCURED') 

125 TYPE 124 

126 DO 128 1=10,70 

TYPE 129,PHI(I,3),PHI(I,5),PHI(I,7) ,PHI(I,9) ,PHI(I, 11) , 

1 PHI(I,13),PHI(I,15),PHI(I,17),PHI(I,19),PHI(I,21) 

128 CONTINUE 

129 FORMAT(10D12.6) 

130 CALL VPLOT 
GO TO 160 

140 OPEN(UNIT=l,FILE='TRAJ.DAT') 

READ (1,4) KK 
DO 150 K=1,KK 
READ(1,2) X0(K),VX0(K) 

READ(1,2) Y0(K),VY0(K) 

READ (1,2) XT(K),YT(K) 
type 155 ,xO(k) ,xt (k) 
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155 formate' x0=',dl2.6 xt=',dl2.6) 

150 CONTINUE 

GO TO 9 
160 STOP 

END 

SUBROUTINE GAUSS(N) 

IMPLICIT REAL *8(A-H,0-Z) 

COMMON/COEFF/VI(5,6) ,AL(5) 

DO 20 K=1,N-1 
B=VI(K,K) 

DO 5 J=K,N+1 
VI(K,J)=VI(K,J)/B 
5 CONTINUE 

DO 15 I=K+1,N 
B=VI(I,K) 

DO 10 J=K,N+1 
V1(I,J)=VI(I,J)-B*VI(K,J) 

10 CONTINUE 

15 CONTINUE 

20 CONTINUE 

AL(N)=VI(N,N+1)/VI(N,N) 

K=N 

25 K=K-1 

AL(K)=VI(K,N+1) 

DO 30 J=K+1,N 
AL(K)=AL(K)-VI(K,J)*AL(J) 

30 CONTINUE 

IF(K.GT.l) GO TO 25 

RETURN 

END 

SUBROUTINE NEOTON ( X , Y , EX , EY , PS I , IK) 

IMPLICIT REAL *8(A-H,0-Z) 

COMMON/COORD/PHI(82,32),XX(82) ,YY(32) ,DX,DY 
DIMENSION Fl-UA.A.Al.FMNCA.^.Al.Xf-ICA) ,:<M1 (A) ,YN(4) ,YN1(A) 
5 FORMAT(2D12.6) 

IF(X.GT.XX(82)) GO TO 220 
IF(X.LT.XX(D) GO TO 220 
IF(Y.GT.YY(32)) GO TO 225 
IF(Y.LT.YYd)) GO TO 220 
HX=DX/80.0 
HY=DY/20.0 
1=0 

10 1 = 1+1 

IF(X.GT.XX(D) GO TO 10 
IF(I.LE.3) GO TO 12 
IB=I-2 
GO TO lA 
12 IB=2 

lA IF(I.EQ.82) GO TO 16 

IE=I+1 
CO TO 18 
IE=I 
J=0 


16 

18 
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J=J+1 

IF(Y.GT.YY(J) ) GO TO 20 
IF(J.LE.3) GO TO 25 
JB=J-2 
GO TO 30 
25 JB=2 

30 IF(J.EQ.32) GO TO 35 

JE=J+1 
GO TO AO 
35 JE=J 

40 1-5M=IE-IB+1 

XM(1)=1.0 
XM1(1)=0.0 
DO 70 M=2,1!M 
IM=IR+M-2 

XM(M)=XM(H-1)*(X-XX(IM)) 

XM1(M)=0.0 

DO 60 1=2, M 

Xl=1.0 

DO 50 K=2,M 

IF(K.EO.I) GO TO 50 

Xl=Xl*(X-XX(IB+K-2)) 

50 CONTINUE 

XM1(M)=XM1(M)+X1 
60 CONTINUE 

70 CONTINUE 

M=JE-JB+1 
YN(l)=1.0 
YN1(1)=0.0 
DO 100 N=2,MN 
JN=JB+N-2 

YN(N)=YN(N-1)*(Y-YY(JN)) 

YN1(N)=0.0 
DO 90 J=2,M 
Yl=1.0 
DO 80 L=2,N 
IF(L.EQ.J) GO TO 80 
Y1=Y 1* (Y-YY ( JB+L-2) ) 

80 CONTINUE 

YN1(N)=YN1(N)+Y1 
90 CONTINUE 

100 CONTINUE 

DO 110 M=1,1-IM 
DO 105 N=1,NN 

FM ( 1 ,M , N ) =PHI ( IB+M- 1 , JB+N- 1 ) 

105 CONTINUE 

110 CONTINUE 

DO 140 L=2,NN 
DO 130 N=L,NN 
DO 120 M=1,MM 

FM(L,H,M) = (FM(L-1,M,N-1)-FM(L-1,M,N))/(HY*(1-D) 
120 CONTINUE 

130 CONTINUE 



140 CONTIMUE 

DO 160 M=1,MM 
DO 150 N=1,NN 
FMM(1,M,N)=FM(N,M,N) 

150 CONTINUE 

160 CONTINUE 

DO 190 K=2,MM 

DO 180 M=K,MM 

DO 170 N=1,NN 

FMN(K,M,N)=(FMN(K-1,M-1,N)-FMN(K-1,M,N))/(HX*(1-K)) 

170 CONTINUE 

180 CONTINUE 

190 CONTINUE 

EX=0.0 
EY=0.0 
PSI=0.0 
DO 210 N=1,NN 

DO 200 M=1,MM 

EX=EX-XM 1 (M ) *YN ( N ) *FMN (M , M , N ) 

EY=EY-YN 1 ( N ) *XH (M) *FMN (M ,M , N ) 
PSI=PSI+XM(M)*YN(N)*FMN(M,M,N) 

200 CONTINUE 

210 CONTINUE 

GO TO 230 
220 IK=-1 

GO TO 230 
225 IK=2 

230 RETURN 

END 

SUBROUTINE DVOGEL(K,KT,KIK) 

IMPLICIT REAL *8(A-H,0-Z) 

COMMON/COORD/PHI(82,32) ,XX(82) ,YY(32) ,DX,DY 
COMMON/DATA/XO(50) ,Y0(50) ,VX0(50) ,VY0(50) ,XT(50) ,YT(50) 
COMHON/SAVE/X(110) , Y( 1 10) ,DT( 1 10) ,KI(50) 
DIM=DSQRT(DX**2+DY**2) 

QM=1.76D11 

IK=0 

IF(KT.GT.l) GO TO 200 
VX=VX0(K) 

VY=VY0(K) 

X(1)=X0(K) 

Y(1)=Y0(K) 

100 CALL NEWTON(X(KT) ,Y(KT) ,EX,EY,PSI,IK) 

101 FORMAK' EX=',D12.4,'EY=',D12.4) 

IF(IK.NE.O) GO TO 400 
AX=QM*EX 

AY=QM*EY 

XH=X ( KT ) - ( 0 . 5*VX-0 . 1 2 5* AX*DT ( KT ) ) *DT ( KT ) 

. YH=Y ( KT ) - ( 0 . 5* VY-0 . 1 25*AY*DT ( KT ) ) *DT ( KT ) 

CALL NEUTON ( XH , Yll , EX , EY , PS I , IK) 

IF(IK.NE.O) GO TO 400 

AXH=OM*EX 

AYH=QM*EY 



200 XH=X(KT)+(0.5*VX+(4.0*AX-AXH)*DT(KT)/24.0)*DT(KT) 

YH = Y ( KT ) + ( 0 . 5 * VY+ ( 4 . 0 * A Y - AYH ) *D T ( KT ) / 2 4 . 0 ) * D T ( KT ) 
CALL NEVrrON ( XH , YH , EX , EY , PSI , IK) 

IF(IK.NE.O) GO TO 400 

AXH=QM*EX 

AYH=QM*EY 

X(KT+1)=X(KT)+(VX+(AX+2.0*AXH)*DT(KT)/6.0)*DT(KT) 
Y (KT+1 ) =Y (KT) +( VY+( AY+2 . 0*AYH) *DT ( KT ) /6 . 0) *DT ( KT ) 
260 FORMAT(2I) 

AAX=AX 

AAY=AY 

CALL NEWTON(X(KT+l) ,Y(KT+1) ,EX,EY,PSI,IK) 

IF(IK.NE.O) GO TO 400 

AX=OM*EX 

AY=QM*EY 

A1SQ=AAX**2+AAY**2 
IF(AlSQ.LT.l.D-6) GO TO 300 

DELD=DSQRT((X(KT+1)-X(KT))**2+(Y(KT+1)-Y(KT) )**2) 
275 FORMATC DELD=' ,012.6/ KT=',13) 

IF(DELD.GT. DIM/40) GO TO 320 
IF(DELD.LT. DIM/120) GO TO 330 
280 VX= VX+DT ( KT ) * ( AX+4 . 0* AXH+AAX ) / 6 . 0 

VY=VY+DT ( KT ) * ( AY+4 . 0* AYH+AAY ) / 6 . 0 
DT(KT+1)=DT(KT) 

GO TO 440 

300 V1SQ=VX**2+VY**2 

DVSq=AlSQ*DT(KT)**2 
IF(VlSQ.GT.DVSq*l.D4) GO TO 280 
1K=1 

GO TO 410 

320 DT(KT)=DT(KT)/1.1 

CO TO 100 

330 DT(KT)=DT(KT)*1.1 

GO TO 100 

400 IF(IK.Eq.-l) TYPE 420 

IF(IK.Eq.2) GO TO 402 

GO TO 440 
402 KI(K)=1 

KIK=1 
TYPE 425 
GO TO 440 
410 TYPE 430 

420 FORMAT(' ELECTRON HAS ESCAPED COORDINATE GRID') 

425 F0RI’1AT(' ELECTRON PAS ESCAPED POTENTIAL REGION') 

430 FORMAT(' PARTICLE VELOCITY AND ACCELERATION BOTH 

1 EqUAL ZERO') 

440 RETURN 

END 

SUBROUTINE VPLOT 
IMPLICIT REAL *8 (A-H.O-Z) 

COMMON/COORD/PHI (82, 32) ,XX(82) ,YY(32) ,DX,DY 
DIMENSION X(250) ,Y(250) 

CALL INITT 
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CALL SELINI 

CALL DWINDO ( SNGL ( XX ( 1 ) ) , SNGL ( XX ( 8 2 ) ) , SNGL ( YY ( 1 ) ) , SNGL ( YY ( 2 2 ) ) ) 
II=IMT ( SNGL ( 780*DY/DX) ) 

CALL TOINDOCO, 780, 0,11) 

CALL MO VE A ( SNGL ( XX ( 2 ) ) , SMGL ( YY ( 2 ) ) ) 

CALL DRAWA ( SNGL ( XX( 8 2 ) ) , SNGL ( YY ( 2 ) ) ) 

PHIMIN=1.0D38 
PHIMAX=-1.0D38 
DO 10 1=2,82 
DO 5 J=2,22 

IF(PKI(I,J) .GT.PHIMIN) GO TO 2 
PHIMIN=PHI(I,J) 

2 IF(PHI(I,J) .LT.PHIMAX) GO TO 5 

PHIMAX=PHI(I,J) 

5 CONTINUE 

10 CONTINUE 

K=0 

12 K=K+1 

V=PHIMIN+K*0.2 

L=0 

DO 20 J=2,22 
DO 15 1=3,82 

IF((PHI(I-1,J)-V)*(PHI(I,J)-V) .GE.0.0) GO TO 15 
L=L+1 

X(L)=XX(I-1)+(XX(I)-XX(I-1))*(V-PHI(I-1,J))/(PHI(I,J)- 
1 PHI(I-1,J)) 

Y(L)=YY(J) 

15 CONTINUE 

20 CONTINUE 

DO 30 1=2,82 
DO 25 J=3,22 

IF((PHI(I,J-1)-V)*(PHI(I,J)-V).GE.0.0) GO TO 25 
L=L+1 

X(L)=XX(I) 

Y(L)=YY(J-1)+(YY(J)-YY(J-1))*(V-PHI(I, J-1))/(PHI(I,J)- 
1 PHI(I,J-D) 

25 CONTINUE 

30 CONTINUE 

LMAX=L 
XMIN=XX(82) 

DO 35 L=1,LMAX 

IF(X(L) .GT.XMIN) GO TO 35 

XMIN=X(L) 

LMIN=L 

35 CONTINUE 

XS=X(1) 

YS=Y(1) 

X(1)=X(LMIN) 

Y(1)=Y(LMIM) 

X(LMIN)=XS 

Y(LMIN)=YS 

CALL HOVEA(SNGL(X(D) ,SNGL(Y(1))) 

LL=1 
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40 LL=LL+1 

DSQ=1.0D38 
DO 45 L=LL,LMAX 

DSQL=(X(L)-X(LL-1))**2+(Y(L)-Y(LL-1))**2 

IF(DSQL.GT.DSQ) GO TO 45 

DSQ=DSQL 

LDMIN=L 

45 CONTINUE 

XS=X(LL) 

YS=Y(LL) 

X(LL)=X(LDMIN) 

Y(LL)=Y(LDMIN) 

X(LDMIN)=XS 

Y(LDHIN)=YS 

CALL DRAUA(SNGL(X(LL)) ,SNGL(Y(LL))) 
IF(LL.EO.LMAX) GO TO 50 
GO TO 40 

50 IF(V+0.2.LE.PHI^^AX) GO TO 12 

CALL FIN1TT(800,500) 

RETURN 

END 



